LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for complex16:

Functions

subroutine zgegs (JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK, INFO)
  ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 
subroutine zgegv (JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
  ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 
subroutine zgees (JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS, LDVS, WORK, LWORK, RWORK, BWORK, INFO)
  ZGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices More...
 
subroutine zgeesx (JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK, INFO)
  ZGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices More...
 
subroutine zgeev (JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
  ZGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 
subroutine zgeevx (BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, RWORK, INFO)
  ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 
subroutine zgges (JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK, BWORK, INFO)
  ZGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices More...
 
subroutine zgges3 (JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK, BWORK, INFO)
  ZGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm) More...
 
subroutine zggesx (JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
  ZGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices More...
 
subroutine zggev (JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
  ZGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 
subroutine zggev3 (JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
  ZGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm) More...
 
subroutine zggevx (BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, BWORK, INFO)
  ZGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 

Detailed Description

This is the group of complex16 eigenvalue driver functions for GE matrices

Function Documentation

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

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

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

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

 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.
          The eigenvalue W(j) is selected if SELECT(W(j)) is true.
[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 has been 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; if
          JOBVS = 'V', LDVS >= N.
[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).
          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]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 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.
Date
November 2011

Definition at line 199 of file zgees.f.

199 *
200 * -- LAPACK driver routine (version 3.4.0) --
201 * -- LAPACK is a software package provided by Univ. of Tennessee, --
202 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203 * November 2011
204 *
205 * .. Scalar Arguments ..
206  CHARACTER jobvs, sort
207  INTEGER info, lda, ldvs, lwork, n, sdim
208 * ..
209 * .. Array Arguments ..
210  LOGICAL bwork( * )
211  DOUBLE PRECISION rwork( * )
212  COMPLEX*16 a( lda, * ), vs( ldvs, * ), w( * ), work( * )
213 * ..
214 * .. Function Arguments ..
215  LOGICAL select
216  EXTERNAL SELECT
217 * ..
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222  DOUBLE PRECISION zero, one
223  parameter( zero = 0.0d0, one = 1.0d0 )
224 * ..
225 * .. Local Scalars ..
226  LOGICAL lquery, scalea, wantst, wantvs
227  INTEGER hswork, i, ibal, icond, ierr, ieval, ihi, ilo,
228  $ itau, iwrk, maxwrk, minwrk
229  DOUBLE PRECISION anrm, bignum, cscale, eps, s, sep, smlnum
230 * ..
231 * .. Local Arrays ..
232  DOUBLE PRECISION dum( 1 )
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL dlabad, xerbla, zcopy, zgebak, zgebal, zgehrd,
237 * ..
238 * .. External Functions ..
239  LOGICAL lsame
240  INTEGER ilaenv
241  DOUBLE PRECISION dlamch, zlange
242  EXTERNAL lsame, ilaenv, dlamch, zlange
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC max, sqrt
246 * ..
247 * .. Executable Statements ..
248 *
249 * Test the input arguments
250 *
251  info = 0
252  lquery = ( lwork.EQ.-1 )
253  wantvs = lsame( jobvs, 'V' )
254  wantst = lsame( sort, 'S' )
255  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
256  info = -1
257  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
258  info = -2
259  ELSE IF( n.LT.0 ) THEN
260  info = -4
261  ELSE IF( lda.LT.max( 1, n ) ) THEN
262  info = -6
263  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
264  info = -10
265  END IF
266 *
267 * Compute workspace
268 * (Note: Comments in the code beginning "Workspace:" describe the
269 * minimal amount of workspace needed at that point in the code,
270 * as well as the preferred amount for good performance.
271 * CWorkspace refers to complex workspace, and RWorkspace to real
272 * workspace. NB refers to the optimal block size for the
273 * immediately following subroutine, as returned by ILAENV.
274 * HSWORK refers to the workspace preferred by ZHSEQR, as
275 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
276 * the worst case.)
277 *
278  IF( info.EQ.0 ) THEN
279  IF( n.EQ.0 ) THEN
280  minwrk = 1
281  maxwrk = 1
282  ELSE
283  maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
284  minwrk = 2*n
285 *
286  CALL zhseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
287  $ work, -1, ieval )
288  hswork = work( 1 )
289 *
290  IF( .NOT.wantvs ) THEN
291  maxwrk = max( maxwrk, hswork )
292  ELSE
293  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
294  $ ' ', n, 1, n, -1 ) )
295  maxwrk = max( maxwrk, hswork )
296  END IF
297  END IF
298  work( 1 ) = maxwrk
299 *
300  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
301  info = -12
302  END IF
303  END IF
304 *
305  IF( info.NE.0 ) THEN
306  CALL xerbla( 'ZGEES ', -info )
307  RETURN
308  ELSE IF( lquery ) THEN
309  RETURN
310  END IF
311 *
312 * Quick return if possible
313 *
314  IF( n.EQ.0 ) THEN
315  sdim = 0
316  RETURN
317  END IF
318 *
319 * Get machine constants
320 *
321  eps = dlamch( 'P' )
322  smlnum = dlamch( 'S' )
323  bignum = one / smlnum
324  CALL dlabad( smlnum, bignum )
325  smlnum = sqrt( smlnum ) / eps
326  bignum = one / smlnum
327 *
328 * Scale A if max element outside range [SMLNUM,BIGNUM]
329 *
330  anrm = zlange( 'M', n, n, a, lda, dum )
331  scalea = .false.
332  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
333  scalea = .true.
334  cscale = smlnum
335  ELSE IF( anrm.GT.bignum ) THEN
336  scalea = .true.
337  cscale = bignum
338  END IF
339  IF( scalea )
340  $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
341 *
342 * Permute the matrix to make it more nearly triangular
343 * (CWorkspace: none)
344 * (RWorkspace: need N)
345 *
346  ibal = 1
347  CALL zgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
348 *
349 * Reduce to upper Hessenberg form
350 * (CWorkspace: need 2*N, prefer N+N*NB)
351 * (RWorkspace: none)
352 *
353  itau = 1
354  iwrk = n + itau
355  CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
356  $ lwork-iwrk+1, ierr )
357 *
358  IF( wantvs ) THEN
359 *
360 * Copy Householder vectors to VS
361 *
362  CALL zlacpy( 'L', n, n, a, lda, vs, ldvs )
363 *
364 * Generate unitary matrix in VS
365 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
366 * (RWorkspace: none)
367 *
368  CALL zunghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
369  $ lwork-iwrk+1, ierr )
370  END IF
371 *
372  sdim = 0
373 *
374 * Perform QR iteration, accumulating Schur vectors in VS if desired
375 * (CWorkspace: need 1, prefer HSWORK (see comments) )
376 * (RWorkspace: none)
377 *
378  iwrk = itau
379  CALL zhseqr( 'S', jobvs, n, ilo, ihi, a, lda, w, vs, ldvs,
380  $ work( iwrk ), lwork-iwrk+1, ieval )
381  IF( ieval.GT.0 )
382  $ info = ieval
383 *
384 * Sort eigenvalues if desired
385 *
386  IF( wantst .AND. info.EQ.0 ) THEN
387  IF( scalea )
388  $ CALL zlascl( 'G', 0, 0, cscale, anrm, n, 1, w, n, ierr )
389  DO 10 i = 1, n
390  bwork( i ) = SELECT( w( i ) )
391  10 CONTINUE
392 *
393 * Reorder eigenvalues and transform Schur vectors
394 * (CWorkspace: none)
395 * (RWorkspace: none)
396 *
397  CALL ztrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, w, sdim,
398  $ s, sep, work( iwrk ), lwork-iwrk+1, icond )
399  END IF
400 *
401  IF( wantvs ) THEN
402 *
403 * Undo balancing
404 * (CWorkspace: none)
405 * (RWorkspace: need N)
406 *
407  CALL zgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs, ldvs,
408  $ ierr )
409  END IF
410 *
411  IF( scalea ) THEN
412 *
413 * Undo scaling for the Schur form of A
414 *
415  CALL zlascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
416  CALL zcopy( n, a, lda+1, w, 1 )
417  END IF
418 *
419  work( 1 ) = maxwrk
420  RETURN
421 *
422 * End of ZGEES
423 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
Definition: zgebal.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:301
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:266
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:128
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:117
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
Definition: zgebak.f:133
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

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 procedure) 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.
Date
November 2011

Definition at line 241 of file zgeesx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeev ( character  JOBVL,
character  JOBVR,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  W,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
 eigenvalues and, optionally, the left and/or right eigenvectors.

 The right eigenvector v(j) of A satisfies
                  A * v(j) = lambda(j) * v(j)
 where lambda(j) is its eigenvalue.
 The left eigenvector u(j) of A satisfies
               u(j)**H * A = lambda(j) * u(j)**H
 where u(j)**H denotes the conjugate transpose of u(j).

 The computed eigenvectors are normalized to have Euclidean norm
 equal to 1 and largest component real.
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N': left eigenvectors of A are not computed;
          = 'V': left eigenvectors of are computed.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N': right eigenvectors of A are not computed;
          = 'V': right eigenvectors of A are computed.
[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 has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]W
          W is COMPLEX*16 array, dimension (N)
          W contains the computed eigenvalues.
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored one
          after another in the columns of VL, in the same order
          as their eigenvalues.
          If JOBVL = 'N', VL is not referenced.
          u(j) = VL(:,j), the j-th column of VL.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.  LDVL >= 1; if
          JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors v(j) are stored one
          after another in the columns of VR, in the same order
          as their eigenvalues.
          If JOBVR = 'N', VR is not referenced.
          v(j) = VR(:,j), the j-th column of VR.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.  LDVR >= 1; if
          JOBVR = 'V', LDVR >= N.
[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).
          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*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, the QR algorithm failed to compute all the
                eigenvalues, and no eigenvectors have been computed;
                elements and i+1:N of W contain eigenvalues which have
                converged.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 179 of file zgeev.f.

179 *
180 * -- LAPACK driver routine (version 3.4.0) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * November 2011
184 *
185 * .. Scalar Arguments ..
186  CHARACTER jobvl, jobvr
187  INTEGER info, lda, ldvl, ldvr, lwork, n
188 * ..
189 * .. Array Arguments ..
190  DOUBLE PRECISION rwork( * )
191  COMPLEX*16 a( lda, * ), vl( ldvl, * ), vr( ldvr, * ),
192  $ w( * ), work( * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  DOUBLE PRECISION zero, one
199  parameter( zero = 0.0d0, one = 1.0d0 )
200 * ..
201 * .. Local Scalars ..
202  LOGICAL lquery, scalea, wantvl, wantvr
203  CHARACTER side
204  INTEGER hswork, i, ibal, ierr, ihi, ilo, irwork, itau,
205  $ iwrk, k, maxwrk, minwrk, nout
206  DOUBLE PRECISION anrm, bignum, cscale, eps, scl, smlnum
207  COMPLEX*16 tmp
208 * ..
209 * .. Local Arrays ..
210  LOGICAL select( 1 )
211  DOUBLE PRECISION dum( 1 )
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL dlabad, xerbla, zdscal, zgebak, zgebal, zgehrd,
216 * ..
217 * .. External Functions ..
218  LOGICAL lsame
219  INTEGER idamax, ilaenv
220  DOUBLE PRECISION dlamch, dznrm2, zlange
221  EXTERNAL lsame, idamax, ilaenv, dlamch, dznrm2, zlange
222 * ..
223 * .. Intrinsic Functions ..
224  INTRINSIC dble, dcmplx, dconjg, dimag, max, sqrt
225 * ..
226 * .. Executable Statements ..
227 *
228 * Test the input arguments
229 *
230  info = 0
231  lquery = ( lwork.EQ.-1 )
232  wantvl = lsame( jobvl, 'V' )
233  wantvr = lsame( jobvr, 'V' )
234  IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
235  info = -1
236  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
237  info = -2
238  ELSE IF( n.LT.0 ) THEN
239  info = -3
240  ELSE IF( lda.LT.max( 1, n ) ) THEN
241  info = -5
242  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
243  info = -8
244  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
245  info = -10
246  END IF
247 *
248 * Compute workspace
249 * (Note: Comments in the code beginning "Workspace:" describe the
250 * minimal amount of workspace needed at that point in the code,
251 * as well as the preferred amount for good performance.
252 * CWorkspace refers to complex workspace, and RWorkspace to real
253 * workspace. NB refers to the optimal block size for the
254 * immediately following subroutine, as returned by ILAENV.
255 * HSWORK refers to the workspace preferred by ZHSEQR, as
256 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
257 * the worst case.)
258 *
259  IF( info.EQ.0 ) THEN
260  IF( n.EQ.0 ) THEN
261  minwrk = 1
262  maxwrk = 1
263  ELSE
264  maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
265  minwrk = 2*n
266  IF( wantvl ) THEN
267  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
268  $ ' ', n, 1, n, -1 ) )
269  CALL zhseqr( 'S', 'V', n, 1, n, a, lda, w, vl, ldvl,
270  $ work, -1, info )
271  ELSE IF( wantvr ) THEN
272  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
273  $ ' ', n, 1, n, -1 ) )
274  CALL zhseqr( 'S', 'V', n, 1, n, a, lda, w, vr, ldvr,
275  $ work, -1, info )
276  ELSE
277  CALL zhseqr( 'E', 'N', n, 1, n, a, lda, w, vr, ldvr,
278  $ work, -1, info )
279  END IF
280  hswork = work( 1 )
281  maxwrk = max( maxwrk, hswork, minwrk )
282  END IF
283  work( 1 ) = maxwrk
284 *
285  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
286  info = -12
287  END IF
288  END IF
289 *
290  IF( info.NE.0 ) THEN
291  CALL xerbla( 'ZGEEV ', -info )
292  RETURN
293  ELSE IF( lquery ) THEN
294  RETURN
295  END IF
296 *
297 * Quick return if possible
298 *
299  IF( n.EQ.0 )
300  $ RETURN
301 *
302 * Get machine constants
303 *
304  eps = dlamch( 'P' )
305  smlnum = dlamch( 'S' )
306  bignum = one / smlnum
307  CALL dlabad( smlnum, bignum )
308  smlnum = sqrt( smlnum ) / eps
309  bignum = one / smlnum
310 *
311 * Scale A if max element outside range [SMLNUM,BIGNUM]
312 *
313  anrm = zlange( 'M', n, n, a, lda, dum )
314  scalea = .false.
315  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
316  scalea = .true.
317  cscale = smlnum
318  ELSE IF( anrm.GT.bignum ) THEN
319  scalea = .true.
320  cscale = bignum
321  END IF
322  IF( scalea )
323  $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
324 *
325 * Balance the matrix
326 * (CWorkspace: none)
327 * (RWorkspace: need N)
328 *
329  ibal = 1
330  CALL zgebal( 'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
331 *
332 * Reduce to upper Hessenberg form
333 * (CWorkspace: need 2*N, prefer N+N*NB)
334 * (RWorkspace: none)
335 *
336  itau = 1
337  iwrk = itau + n
338  CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
339  $ lwork-iwrk+1, ierr )
340 *
341  IF( wantvl ) THEN
342 *
343 * Want left eigenvectors
344 * Copy Householder vectors to VL
345 *
346  side = 'L'
347  CALL zlacpy( 'L', n, n, a, lda, vl, ldvl )
348 *
349 * Generate unitary matrix in VL
350 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
351 * (RWorkspace: none)
352 *
353  CALL zunghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
354  $ lwork-iwrk+1, ierr )
355 *
356 * Perform QR iteration, accumulating Schur vectors in VL
357 * (CWorkspace: need 1, prefer HSWORK (see comments) )
358 * (RWorkspace: none)
359 *
360  iwrk = itau
361  CALL zhseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vl, ldvl,
362  $ work( iwrk ), lwork-iwrk+1, info )
363 *
364  IF( wantvr ) THEN
365 *
366 * Want left and right eigenvectors
367 * Copy Schur vectors to VR
368 *
369  side = 'B'
370  CALL zlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
371  END IF
372 *
373  ELSE IF( wantvr ) THEN
374 *
375 * Want right eigenvectors
376 * Copy Householder vectors to VR
377 *
378  side = 'R'
379  CALL zlacpy( 'L', n, n, a, lda, vr, ldvr )
380 *
381 * Generate unitary matrix in VR
382 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
383 * (RWorkspace: none)
384 *
385  CALL zunghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
386  $ lwork-iwrk+1, ierr )
387 *
388 * Perform QR iteration, accumulating Schur vectors in VR
389 * (CWorkspace: need 1, prefer HSWORK (see comments) )
390 * (RWorkspace: none)
391 *
392  iwrk = itau
393  CALL zhseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vr, ldvr,
394  $ work( iwrk ), lwork-iwrk+1, info )
395 *
396  ELSE
397 *
398 * Compute eigenvalues only
399 * (CWorkspace: need 1, prefer HSWORK (see comments) )
400 * (RWorkspace: none)
401 *
402  iwrk = itau
403  CALL zhseqr( 'E', 'N', n, ilo, ihi, a, lda, w, vr, ldvr,
404  $ work( iwrk ), lwork-iwrk+1, info )
405  END IF
406 *
407 * If INFO > 0 from ZHSEQR, then quit
408 *
409  IF( info.GT.0 )
410  $ GO TO 50
411 *
412  IF( wantvl .OR. wantvr ) THEN
413 *
414 * Compute left and/or right eigenvectors
415 * (CWorkspace: need 2*N)
416 * (RWorkspace: need 2*N)
417 *
418  irwork = ibal + n
419  CALL ztrevc( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
420  $ n, nout, work( iwrk ), rwork( irwork ), ierr )
421  END IF
422 *
423  IF( wantvl ) THEN
424 *
425 * Undo balancing of left eigenvectors
426 * (CWorkspace: none)
427 * (RWorkspace: need N)
428 *
429  CALL zgebak( 'B', 'L', n, ilo, ihi, rwork( ibal ), n, vl, ldvl,
430  $ ierr )
431 *
432 * Normalize left eigenvectors and make largest component real
433 *
434  DO 20 i = 1, n
435  scl = one / dznrm2( n, vl( 1, i ), 1 )
436  CALL zdscal( n, scl, vl( 1, i ), 1 )
437  DO 10 k = 1, n
438  rwork( irwork+k-1 ) = dble( vl( k, i ) )**2 +
439  $ dimag( vl( k, i ) )**2
440  10 CONTINUE
441  k = idamax( n, rwork( irwork ), 1 )
442  tmp = dconjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
443  CALL zscal( n, tmp, vl( 1, i ), 1 )
444  vl( k, i ) = dcmplx( dble( vl( k, i ) ), zero )
445  20 CONTINUE
446  END IF
447 *
448  IF( wantvr ) THEN
449 *
450 * Undo balancing of right eigenvectors
451 * (CWorkspace: none)
452 * (RWorkspace: need N)
453 *
454  CALL zgebak( 'B', 'R', n, ilo, ihi, rwork( ibal ), n, vr, ldvr,
455  $ ierr )
456 *
457 * Normalize right eigenvectors and make largest component real
458 *
459  DO 40 i = 1, n
460  scl = one / dznrm2( n, vr( 1, i ), 1 )
461  CALL zdscal( n, scl, vr( 1, i ), 1 )
462  DO 30 k = 1, n
463  rwork( irwork+k-1 ) = dble( vr( k, i ) )**2 +
464  $ dimag( vr( k, i ) )**2
465  30 CONTINUE
466  k = idamax( n, rwork( irwork ), 1 )
467  tmp = dconjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
468  CALL zscal( n, tmp, vr( 1, i ), 1 )
469  vr( k, i ) = dcmplx( dble( vr( k, i ) ), zero )
470  40 CONTINUE
471  END IF
472 *
473 * Undo scaling if necessary
474 *
475  50 CONTINUE
476  IF( scalea ) THEN
477  CALL zlascl( 'G', 0, 0, cscale, anrm, n-info, 1, w( info+1 ),
478  $ max( n-info, 1 ), ierr )
479  IF( info.GT.0 ) THEN
480  CALL zlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, w, n, ierr )
481  END IF
482  END IF
483 *
484  work( 1 ) = maxwrk
485  RETURN
486 *
487 * End of ZGEEV
488 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
Definition: zgebal.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine ztrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTREVC
Definition: ztrevc.f:220
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:301
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:128
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:117
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
Definition: zgebak.f:133
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeevx ( character  BALANC,
character  JOBVL,
character  JOBVR,
character  SENSE,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  W,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
integer  ILO,
integer  IHI,
double precision, dimension( * )  SCALE,
double precision  ABNRM,
double precision, dimension( * )  RCONDE,
double precision, dimension( * )  RCONDV,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 ZGEEVX computes for an N-by-N complex nonsymmetric matrix A, the
 eigenvalues and, optionally, the left and/or right eigenvectors.

 Optionally also, it computes a balancing transformation to improve
 the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
 SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
 (RCONDE), and reciprocal condition numbers for the right
 eigenvectors (RCONDV).

 The right eigenvector v(j) of A satisfies
                  A * v(j) = lambda(j) * v(j)
 where lambda(j) is its eigenvalue.
 The left eigenvector u(j) of A satisfies
               u(j)**H * A = lambda(j) * u(j)**H
 where u(j)**H denotes the conjugate transpose of u(j).

 The computed eigenvectors are normalized to have Euclidean norm
 equal to 1 and largest component real.

 Balancing a matrix means permuting the rows and columns to make it
 more nearly upper triangular, and applying a diagonal similarity
 transformation D * A * D**(-1), where D is a diagonal matrix, to
 make its rows and columns closer in norm and the condition numbers
 of its eigenvalues and eigenvectors smaller.  The computed
 reciprocal condition numbers correspond to the balanced matrix.
 Permuting rows and columns will not change the condition numbers
 (in exact arithmetic) but diagonal scaling will.  For further
 explanation of balancing, see section 4.10.2 of the LAPACK
 Users' Guide.
Parameters
[in]BALANC
          BALANC is CHARACTER*1
          Indicates how the input matrix should be diagonally scaled
          and/or permuted to improve the conditioning of its
          eigenvalues.
          = 'N': Do not diagonally scale or permute;
          = 'P': Perform permutations to make the matrix more nearly
                 upper triangular. Do not diagonally scale;
          = 'S': Diagonally scale the matrix, ie. replace A by
                 D*A*D**(-1), where D is a diagonal matrix chosen
                 to make the rows and columns of A more equal in
                 norm. Do not permute;
          = 'B': Both diagonally scale and permute A.

          Computed reciprocal condition numbers will be for the matrix
          after balancing and/or permuting. Permuting does not change
          condition numbers (in exact arithmetic), but balancing does.
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N': left eigenvectors of A are not computed;
          = 'V': left eigenvectors of A are computed.
          If SENSE = 'E' or 'B', JOBVL must = 'V'.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N': right eigenvectors of A are not computed;
          = 'V': right eigenvectors of A are computed.
          If SENSE = 'E' or 'B', JOBVR must = 'V'.
[in]SENSE
          SENSE is CHARACTER*1
          Determines which reciprocal condition numbers are computed.
          = 'N': None are computed;
          = 'E': Computed for eigenvalues only;
          = 'V': Computed for right eigenvectors only;
          = 'B': Computed for eigenvalues and right eigenvectors.

          If SENSE = 'E' or 'B', both left and right eigenvectors
          must also be computed (JOBVL = 'V' and JOBVR = 'V').
[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 has been overwritten.  If JOBVL = 'V' or
          JOBVR = 'V', A contains the Schur form of the balanced
          version of the matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]W
          W is COMPLEX*16 array, dimension (N)
          W contains the computed eigenvalues.
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored one
          after another in the columns of VL, in the same order
          as their eigenvalues.
          If JOBVL = 'N', VL is not referenced.
          u(j) = VL(:,j), the j-th column of VL.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.  LDVL >= 1; if
          JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors v(j) are stored one
          after another in the columns of VR, in the same order
          as their eigenvalues.
          If JOBVR = 'N', VR is not referenced.
          v(j) = VR(:,j), the j-th column of VR.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.  LDVR >= 1; if
          JOBVR = 'V', LDVR >= N.
[out]ILO
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are integer values determined when A was
          balanced.  The balanced A(i,j) = 0 if I > J and
          J = 1,...,ILO-1 or I = IHI+1,...,N.
[out]SCALE
          SCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          when balancing A.  If P(j) is the index of the row and column
          interchanged with row and column j, and D(j) is the scaling
          factor applied to row and column j, then
          SCALE(J) = P(J),    for J = 1,...,ILO-1
                   = D(J),    for J = ILO,...,IHI
                   = P(J)     for J = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]ABNRM
          ABNRM is DOUBLE PRECISION
          The one-norm of the balanced matrix (the maximum
          of the sum of absolute values of elements of any column).
[out]RCONDE
          RCONDE is DOUBLE PRECISION array, dimension (N)
          RCONDE(j) is the reciprocal condition number of the j-th
          eigenvalue.
[out]RCONDV
          RCONDV is DOUBLE PRECISION array, dimension (N)
          RCONDV(j) is the reciprocal condition number of the j-th
          right eigenvector.
[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.  If SENSE = 'N' or 'E',
          LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
          LWORK >= N*N+2*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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*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, the QR algorithm failed to compute all the
                eigenvalues, and no eigenvectors or condition numbers
                have been computed; elements 1:ILO-1 and i+1:N of W
                contain eigenvalues which have converged.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 287 of file zgeevx.f.

287 *
288 * -- LAPACK driver routine (version 3.4.0) --
289 * -- LAPACK is a software package provided by Univ. of Tennessee, --
290 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
291 * November 2011
292 *
293 * .. Scalar Arguments ..
294  CHARACTER balanc, jobvl, jobvr, sense
295  INTEGER ihi, ilo, info, lda, ldvl, ldvr, lwork, n
296  DOUBLE PRECISION abnrm
297 * ..
298 * .. Array Arguments ..
299  DOUBLE PRECISION rconde( * ), rcondv( * ), rwork( * ),
300  $ scale( * )
301  COMPLEX*16 a( lda, * ), vl( ldvl, * ), vr( ldvr, * ),
302  $ w( * ), work( * )
303 * ..
304 *
305 * =====================================================================
306 *
307 * .. Parameters ..
308  DOUBLE PRECISION zero, one
309  parameter( zero = 0.0d0, one = 1.0d0 )
310 * ..
311 * .. Local Scalars ..
312  LOGICAL lquery, scalea, wantvl, wantvr, wntsnb, wntsne,
313  $ wntsnn, wntsnv
314  CHARACTER job, side
315  INTEGER hswork, i, icond, ierr, itau, iwrk, k, maxwrk,
316  $ minwrk, nout
317  DOUBLE PRECISION anrm, bignum, cscale, eps, scl, smlnum
318  COMPLEX*16 tmp
319 * ..
320 * .. Local Arrays ..
321  LOGICAL select( 1 )
322  DOUBLE PRECISION dum( 1 )
323 * ..
324 * .. External Subroutines ..
325  EXTERNAL dlabad, dlascl, xerbla, zdscal, zgebak, zgebal,
327  $ ztrsna, zunghr
328 * ..
329 * .. External Functions ..
330  LOGICAL lsame
331  INTEGER idamax, ilaenv
332  DOUBLE PRECISION dlamch, dznrm2, zlange
333  EXTERNAL lsame, idamax, ilaenv, dlamch, dznrm2, zlange
334 * ..
335 * .. Intrinsic Functions ..
336  INTRINSIC dble, dcmplx, dconjg, dimag, max, sqrt
337 * ..
338 * .. Executable Statements ..
339 *
340 * Test the input arguments
341 *
342  info = 0
343  lquery = ( lwork.EQ.-1 )
344  wantvl = lsame( jobvl, 'V' )
345  wantvr = lsame( jobvr, 'V' )
346  wntsnn = lsame( sense, 'N' )
347  wntsne = lsame( sense, 'E' )
348  wntsnv = lsame( sense, 'V' )
349  wntsnb = lsame( sense, 'B' )
350  IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc, 'S' ) .OR.
351  $ lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) ) THEN
352  info = -1
353  ELSE IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
354  info = -2
355  ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
356  info = -3
357  ELSE IF( .NOT.( wntsnn .OR. wntsne .OR. wntsnb .OR. wntsnv ) .OR.
358  $ ( ( wntsne .OR. wntsnb ) .AND. .NOT.( wantvl .AND.
359  $ wantvr ) ) ) THEN
360  info = -4
361  ELSE IF( n.LT.0 ) THEN
362  info = -5
363  ELSE IF( lda.LT.max( 1, n ) ) THEN
364  info = -7
365  ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
366  info = -10
367  ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
368  info = -12
369  END IF
370 *
371 * Compute workspace
372 * (Note: Comments in the code beginning "Workspace:" describe the
373 * minimal amount of workspace needed at that point in the code,
374 * as well as the preferred amount for good performance.
375 * CWorkspace refers to complex workspace, and RWorkspace to real
376 * workspace. NB refers to the optimal block size for the
377 * immediately following subroutine, as returned by ILAENV.
378 * HSWORK refers to the workspace preferred by ZHSEQR, as
379 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
380 * the worst case.)
381 *
382  IF( info.EQ.0 ) THEN
383  IF( n.EQ.0 ) THEN
384  minwrk = 1
385  maxwrk = 1
386  ELSE
387  maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
388 *
389  IF( wantvl ) THEN
390  CALL zhseqr( 'S', 'V', n, 1, n, a, lda, w, vl, ldvl,
391  $ work, -1, info )
392  ELSE IF( wantvr ) THEN
393  CALL zhseqr( 'S', 'V', n, 1, n, a, lda, w, vr, ldvr,
394  $ work, -1, info )
395  ELSE
396  IF( wntsnn ) THEN
397  CALL zhseqr( 'E', 'N', n, 1, n, a, lda, w, vr, ldvr,
398  $ work, -1, info )
399  ELSE
400  CALL zhseqr( 'S', 'N', n, 1, n, a, lda, w, vr, ldvr,
401  $ work, -1, info )
402  END IF
403  END IF
404  hswork = work( 1 )
405 *
406  IF( ( .NOT.wantvl ) .AND. ( .NOT.wantvr ) ) THEN
407  minwrk = 2*n
408  IF( .NOT.( wntsnn .OR. wntsne ) )
409  $ minwrk = max( minwrk, n*n + 2*n )
410  maxwrk = max( maxwrk, hswork )
411  IF( .NOT.( wntsnn .OR. wntsne ) )
412  $ maxwrk = max( maxwrk, n*n + 2*n )
413  ELSE
414  minwrk = 2*n
415  IF( .NOT.( wntsnn .OR. wntsne ) )
416  $ minwrk = max( minwrk, n*n + 2*n )
417  maxwrk = max( maxwrk, hswork )
418  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
419  $ ' ', n, 1, n, -1 ) )
420  IF( .NOT.( wntsnn .OR. wntsne ) )
421  $ maxwrk = max( maxwrk, n*n + 2*n )
422  maxwrk = max( maxwrk, 2*n )
423  END IF
424  maxwrk = max( maxwrk, minwrk )
425  END IF
426  work( 1 ) = maxwrk
427 *
428  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
429  info = -20
430  END IF
431  END IF
432 *
433  IF( info.NE.0 ) THEN
434  CALL xerbla( 'ZGEEVX', -info )
435  RETURN
436  ELSE IF( lquery ) THEN
437  RETURN
438  END IF
439 *
440 * Quick return if possible
441 *
442  IF( n.EQ.0 )
443  $ RETURN
444 *
445 * Get machine constants
446 *
447  eps = dlamch( 'P' )
448  smlnum = dlamch( 'S' )
449  bignum = one / smlnum
450  CALL dlabad( smlnum, bignum )
451  smlnum = sqrt( smlnum ) / eps
452  bignum = one / smlnum
453 *
454 * Scale A if max element outside range [SMLNUM,BIGNUM]
455 *
456  icond = 0
457  anrm = zlange( 'M', n, n, a, lda, dum )
458  scalea = .false.
459  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
460  scalea = .true.
461  cscale = smlnum
462  ELSE IF( anrm.GT.bignum ) THEN
463  scalea = .true.
464  cscale = bignum
465  END IF
466  IF( scalea )
467  $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
468 *
469 * Balance the matrix and compute ABNRM
470 *
471  CALL zgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
472  abnrm = zlange( '1', n, n, a, lda, dum )
473  IF( scalea ) THEN
474  dum( 1 ) = abnrm
475  CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
476  abnrm = dum( 1 )
477  END IF
478 *
479 * Reduce to upper Hessenberg form
480 * (CWorkspace: need 2*N, prefer N+N*NB)
481 * (RWorkspace: none)
482 *
483  itau = 1
484  iwrk = itau + n
485  CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
486  $ lwork-iwrk+1, ierr )
487 *
488  IF( wantvl ) THEN
489 *
490 * Want left eigenvectors
491 * Copy Householder vectors to VL
492 *
493  side = 'L'
494  CALL zlacpy( 'L', n, n, a, lda, vl, ldvl )
495 *
496 * Generate unitary matrix in VL
497 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
498 * (RWorkspace: none)
499 *
500  CALL zunghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
501  $ lwork-iwrk+1, ierr )
502 *
503 * Perform QR iteration, accumulating Schur vectors in VL
504 * (CWorkspace: need 1, prefer HSWORK (see comments) )
505 * (RWorkspace: none)
506 *
507  iwrk = itau
508  CALL zhseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vl, ldvl,
509  $ work( iwrk ), lwork-iwrk+1, info )
510 *
511  IF( wantvr ) THEN
512 *
513 * Want left and right eigenvectors
514 * Copy Schur vectors to VR
515 *
516  side = 'B'
517  CALL zlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
518  END IF
519 *
520  ELSE IF( wantvr ) THEN
521 *
522 * Want right eigenvectors
523 * Copy Householder vectors to VR
524 *
525  side = 'R'
526  CALL zlacpy( 'L', n, n, a, lda, vr, ldvr )
527 *
528 * Generate unitary matrix in VR
529 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
530 * (RWorkspace: none)
531 *
532  CALL zunghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
533  $ lwork-iwrk+1, ierr )
534 *
535 * Perform QR iteration, accumulating Schur vectors in VR
536 * (CWorkspace: need 1, prefer HSWORK (see comments) )
537 * (RWorkspace: none)
538 *
539  iwrk = itau
540  CALL zhseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vr, ldvr,
541  $ work( iwrk ), lwork-iwrk+1, info )
542 *
543  ELSE
544 *
545 * Compute eigenvalues only
546 * If condition numbers desired, compute Schur form
547 *
548  IF( wntsnn ) THEN
549  job = 'E'
550  ELSE
551  job = 'S'
552  END IF
553 *
554 * (CWorkspace: need 1, prefer HSWORK (see comments) )
555 * (RWorkspace: none)
556 *
557  iwrk = itau
558  CALL zhseqr( job, 'N', n, ilo, ihi, a, lda, w, vr, ldvr,
559  $ work( iwrk ), lwork-iwrk+1, info )
560  END IF
561 *
562 * If INFO > 0 from ZHSEQR, then quit
563 *
564  IF( info.GT.0 )
565  $ GO TO 50
566 *
567  IF( wantvl .OR. wantvr ) THEN
568 *
569 * Compute left and/or right eigenvectors
570 * (CWorkspace: need 2*N)
571 * (RWorkspace: need N)
572 *
573  CALL ztrevc( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
574  $ n, nout, work( iwrk ), rwork, ierr )
575  END IF
576 *
577 * Compute condition numbers if desired
578 * (CWorkspace: need N*N+2*N unless SENSE = 'E')
579 * (RWorkspace: need 2*N unless SENSE = 'E')
580 *
581  IF( .NOT.wntsnn ) THEN
582  CALL ztrsna( sense, 'A', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
583  $ rconde, rcondv, n, nout, work( iwrk ), n, rwork,
584  $ icond )
585  END IF
586 *
587  IF( wantvl ) THEN
588 *
589 * Undo balancing of left eigenvectors
590 *
591  CALL zgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,
592  $ ierr )
593 *
594 * Normalize left eigenvectors and make largest component real
595 *
596  DO 20 i = 1, n
597  scl = one / dznrm2( n, vl( 1, i ), 1 )
598  CALL zdscal( n, scl, vl( 1, i ), 1 )
599  DO 10 k = 1, n
600  rwork( k ) = dble( vl( k, i ) )**2 +
601  $ dimag( vl( k, i ) )**2
602  10 CONTINUE
603  k = idamax( n, rwork, 1 )
604  tmp = dconjg( vl( k, i ) ) / sqrt( rwork( k ) )
605  CALL zscal( n, tmp, vl( 1, i ), 1 )
606  vl( k, i ) = dcmplx( dble( vl( k, i ) ), zero )
607  20 CONTINUE
608  END IF
609 *
610  IF( wantvr ) THEN
611 *
612 * Undo balancing of right eigenvectors
613 *
614  CALL zgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,
615  $ ierr )
616 *
617 * Normalize right eigenvectors and make largest component real
618 *
619  DO 40 i = 1, n
620  scl = one / dznrm2( n, vr( 1, i ), 1 )
621  CALL zdscal( n, scl, vr( 1, i ), 1 )
622  DO 30 k = 1, n
623  rwork( k ) = dble( vr( k, i ) )**2 +
624  $ dimag( vr( k, i ) )**2
625  30 CONTINUE
626  k = idamax( n, rwork, 1 )
627  tmp = dconjg( vr( k, i ) ) / sqrt( rwork( k ) )
628  CALL zscal( n, tmp, vr( 1, i ), 1 )
629  vr( k, i ) = dcmplx( dble( vr( k, i ) ), zero )
630  40 CONTINUE
631  END IF
632 *
633 * Undo scaling if necessary
634 *
635  50 CONTINUE
636  IF( scalea ) THEN
637  CALL zlascl( 'G', 0, 0, cscale, anrm, n-info, 1, w( info+1 ),
638  $ max( n-info, 1 ), ierr )
639  IF( info.EQ.0 ) THEN
640  IF( ( wntsnv .OR. wntsnb ) .AND. icond.EQ.0 )
641  $ CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, rcondv, n,
642  $ ierr )
643  ELSE
644  CALL zlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, w, n, ierr )
645  END IF
646  END IF
647 *
648  work( 1 ) = maxwrk
649  RETURN
650 *
651 * End of ZGEEVX
652 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
Definition: zgebal.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine ztrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTREVC
Definition: ztrevc.f:220
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:301
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:128
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine ztrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, INFO)
ZTRSNA
Definition: ztrsna.f:251
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:117
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:141
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
Definition: zgebak.f:133
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgegs ( character  JOBVSL,
character  JOBVSR,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvsl, * )  VSL,
integer  LDVSL,
complex*16, dimension( ldvsr, * )  VSR,
integer  LDVSR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 This routine is deprecated and has been replaced by routine ZGGES.

 ZGEGS computes the eigenvalues, Schur form, and, optionally, the
 left and or/right Schur vectors of a complex matrix pair (A,B).
 Given two square matrices A and B, the generalized Schur
 factorization has the form
 
    A = Q*S*Z**H,  B = Q*T*Z**H
 
 where Q and Z are unitary matrices and S and T are upper triangular.
 The columns of Q are the left Schur vectors
 and the columns of Z are the right Schur vectors.
 
 If only the eigenvalues of (A,B) are needed, the driver routine
 ZGEGV should be used instead.  See ZGEGV for a description of the
 eigenvalues of the generalized nonsymmetric eigenvalue problem
 (GNEP).
Parameters
[in]JOBVSL
          JOBVSL is CHARACTER*1
          = 'N':  do not compute the left Schur vectors;
          = 'V':  compute the left Schur vectors (returned in VSL).
[in]JOBVSR
          JOBVSR is CHARACTER*1
          = 'N':  do not compute the right Schur vectors;
          = 'V':  compute the right Schur vectors (returned in VSR).
[in]N
          N is INTEGER
          The order of the matrices A, B, VSL, and VSR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the matrix A.
          On exit, the upper triangular matrix S from the generalized
          Schur factorization.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the matrix B.
          On exit, the upper triangular matrix T from the generalized
          Schur factorization.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
          The complex scalars alpha that define the eigenvalues of
          GNEP.  ALPHA(j) = S(j,j), the diagonal element of the Schur
          form of A.
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          The non-negative real scalars beta that define the
          eigenvalues of GNEP.  BETA(j) = T(j,j), the diagonal element
          of the triangular factor T.

          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
          represent the j-th eigenvalue of the matrix pair (A,B), in
          one of the forms lambda = alpha/beta or mu = beta/alpha.
          Since either lambda or mu may overflow, they should not,
          in general, be computed.
[out]VSL
          VSL is COMPLEX*16 array, dimension (LDVSL,N)
          If JOBVSL = 'V', the matrix of left Schur vectors Q.
          Not referenced if JOBVSL = 'N'.
[in]LDVSL
          LDVSL is INTEGER
          The leading dimension of the matrix VSL. LDVSL >= 1, and
          if JOBVSL = 'V', LDVSL >= N.
[out]VSR
          VSR is COMPLEX*16 array, dimension (LDVSR,N)
          If JOBVSR = 'V', the matrix of right Schur vectors Z.
          Not referenced if JOBVSR = 'N'.
[in]LDVSR
          LDVSR is INTEGER
          The leading dimension of the matrix VSR. LDVSR >= 1, and
          if JOBVSR = 'V', LDVSR >= N.
[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).
          For good performance, LWORK must generally be larger.
          To compute the optimal value of LWORK, call ILAENV to get
          blocksizes (for ZGEQRF, ZUNMQR, and CUNGQR.)  Then compute:
          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and CUNGQR;
          the optimal LWORK is N*(NB+1).

          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (3*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          =1,...,N:
                The QZ iteration failed.  (A,B) are not in Schur
                form, but ALPHA(j) and BETA(j) should be correct for
                j=INFO+1,...,N.
          > N:  errors that usually indicate LAPACK problems:
                =N+1: error return from ZGGBAL
                =N+2: error return from ZGEQRF
                =N+3: error return from ZUNMQR
                =N+4: error return from ZUNGQR
                =N+5: error return from ZGGHRD
                =N+6: error return from ZHGEQZ (other than failed
                                               iteration)
                =N+7: error return from ZGGBAK (computing VSL)
                =N+8: error return from ZGGBAK (computing VSR)
                =N+9: error return from ZLASCL (various places)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 227 of file zgegs.f.

227 *
228 * -- LAPACK driver routine (version 3.4.0) --
229 * -- LAPACK is a software package provided by Univ. of Tennessee, --
230 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231 * November 2011
232 *
233 * .. Scalar Arguments ..
234  CHARACTER jobvsl, jobvsr
235  INTEGER info, lda, ldb, ldvsl, ldvsr, lwork, n
236 * ..
237 * .. Array Arguments ..
238  DOUBLE PRECISION rwork( * )
239  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
240  $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
241  $ work( * )
242 * ..
243 *
244 * =====================================================================
245 *
246 * .. Parameters ..
247  DOUBLE PRECISION zero, one
248  parameter( zero = 0.0d0, one = 1.0d0 )
249  COMPLEX*16 czero, cone
250  parameter( czero = ( 0.0d0, 0.0d0 ),
251  $ cone = ( 1.0d0, 0.0d0 ) )
252 * ..
253 * .. Local Scalars ..
254  LOGICAL ilascl, ilbscl, ilvsl, ilvsr, lquery
255  INTEGER icols, ihi, iinfo, ijobvl, ijobvr, ileft, ilo,
256  $ iright, irows, irwork, itau, iwork, lopt,
257  $ lwkmin, lwkopt, nb, nb1, nb2, nb3
258  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
259  $ safmin, smlnum
260 * ..
261 * .. External Subroutines ..
262  EXTERNAL xerbla, zgeqrf, zggbak, zggbal, zgghrd, zhgeqz,
264 * ..
265 * .. External Functions ..
266  LOGICAL lsame
267  INTEGER ilaenv
268  DOUBLE PRECISION dlamch, zlange
269  EXTERNAL lsame, ilaenv, dlamch, zlange
270 * ..
271 * .. Intrinsic Functions ..
272  INTRINSIC int, max
273 * ..
274 * .. Executable Statements ..
275 *
276 * Decode the input arguments
277 *
278  IF( lsame( jobvsl, 'N' ) ) THEN
279  ijobvl = 1
280  ilvsl = .false.
281  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
282  ijobvl = 2
283  ilvsl = .true.
284  ELSE
285  ijobvl = -1
286  ilvsl = .false.
287  END IF
288 *
289  IF( lsame( jobvsr, 'N' ) ) THEN
290  ijobvr = 1
291  ilvsr = .false.
292  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
293  ijobvr = 2
294  ilvsr = .true.
295  ELSE
296  ijobvr = -1
297  ilvsr = .false.
298  END IF
299 *
300 * Test the input arguments
301 *
302  lwkmin = max( 2*n, 1 )
303  lwkopt = lwkmin
304  work( 1 ) = lwkopt
305  lquery = ( lwork.EQ.-1 )
306  info = 0
307  IF( ijobvl.LE.0 ) THEN
308  info = -1
309  ELSE IF( ijobvr.LE.0 ) THEN
310  info = -2
311  ELSE IF( n.LT.0 ) THEN
312  info = -3
313  ELSE IF( lda.LT.max( 1, n ) ) THEN
314  info = -5
315  ELSE IF( ldb.LT.max( 1, n ) ) THEN
316  info = -7
317  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
318  info = -11
319  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
320  info = -13
321  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
322  info = -15
323  END IF
324 *
325  IF( info.EQ.0 ) THEN
326  nb1 = ilaenv( 1, 'ZGEQRF', ' ', n, n, -1, -1 )
327  nb2 = ilaenv( 1, 'ZUNMQR', ' ', n, n, n, -1 )
328  nb3 = ilaenv( 1, 'ZUNGQR', ' ', n, n, n, -1 )
329  nb = max( nb1, nb2, nb3 )
330  lopt = n*( nb+1 )
331  work( 1 ) = lopt
332  END IF
333 *
334  IF( info.NE.0 ) THEN
335  CALL xerbla( 'ZGEGS ', -info )
336  RETURN
337  ELSE IF( lquery ) THEN
338  RETURN
339  END IF
340 *
341 * Quick return if possible
342 *
343  IF( n.EQ.0 )
344  $ RETURN
345 *
346 * Get machine constants
347 *
348  eps = dlamch( 'E' )*dlamch( 'B' )
349  safmin = dlamch( 'S' )
350  smlnum = n*safmin / eps
351  bignum = one / smlnum
352 *
353 * Scale A if max element outside range [SMLNUM,BIGNUM]
354 *
355  anrm = zlange( 'M', n, n, a, lda, rwork )
356  ilascl = .false.
357  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
358  anrmto = smlnum
359  ilascl = .true.
360  ELSE IF( anrm.GT.bignum ) THEN
361  anrmto = bignum
362  ilascl = .true.
363  END IF
364 *
365  IF( ilascl ) THEN
366  CALL zlascl( 'G', -1, -1, anrm, anrmto, n, n, a, lda, iinfo )
367  IF( iinfo.NE.0 ) THEN
368  info = n + 9
369  RETURN
370  END IF
371  END IF
372 *
373 * Scale B if max element outside range [SMLNUM,BIGNUM]
374 *
375  bnrm = zlange( 'M', n, n, b, ldb, rwork )
376  ilbscl = .false.
377  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
378  bnrmto = smlnum
379  ilbscl = .true.
380  ELSE IF( bnrm.GT.bignum ) THEN
381  bnrmto = bignum
382  ilbscl = .true.
383  END IF
384 *
385  IF( ilbscl ) THEN
386  CALL zlascl( 'G', -1, -1, bnrm, bnrmto, n, n, b, ldb, iinfo )
387  IF( iinfo.NE.0 ) THEN
388  info = n + 9
389  RETURN
390  END IF
391  END IF
392 *
393 * Permute the matrix to make it more nearly triangular
394 *
395  ileft = 1
396  iright = n + 1
397  irwork = iright + n
398  iwork = 1
399  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
400  $ rwork( iright ), rwork( irwork ), iinfo )
401  IF( iinfo.NE.0 ) THEN
402  info = n + 1
403  GO TO 10
404  END IF
405 *
406 * Reduce B to triangular form, and initialize VSL and/or VSR
407 *
408  irows = ihi + 1 - ilo
409  icols = n + 1 - ilo
410  itau = iwork
411  iwork = itau + irows
412  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
413  $ work( iwork ), lwork+1-iwork, iinfo )
414  IF( iinfo.GE.0 )
415  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
416  IF( iinfo.NE.0 ) THEN
417  info = n + 2
418  GO TO 10
419  END IF
420 *
421  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
422  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
423  $ lwork+1-iwork, iinfo )
424  IF( iinfo.GE.0 )
425  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
426  IF( iinfo.NE.0 ) THEN
427  info = n + 3
428  GO TO 10
429  END IF
430 *
431  IF( ilvsl ) THEN
432  CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
433  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
434  $ vsl( ilo+1, ilo ), ldvsl )
435  CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
436  $ work( itau ), work( iwork ), lwork+1-iwork,
437  $ iinfo )
438  IF( iinfo.GE.0 )
439  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
440  IF( iinfo.NE.0 ) THEN
441  info = n + 4
442  GO TO 10
443  END IF
444  END IF
445 *
446  IF( ilvsr )
447  $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
448 *
449 * Reduce to generalized Hessenberg form
450 *
451  CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
452  $ ldvsl, vsr, ldvsr, iinfo )
453  IF( iinfo.NE.0 ) THEN
454  info = n + 5
455  GO TO 10
456  END IF
457 *
458 * Perform QZ algorithm, computing Schur vectors if desired
459 *
460  iwork = itau
461  CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
462  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwork ),
463  $ lwork+1-iwork, rwork( irwork ), iinfo )
464  IF( iinfo.GE.0 )
465  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
466  IF( iinfo.NE.0 ) THEN
467  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
468  info = iinfo
469  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
470  info = iinfo - n
471  ELSE
472  info = n + 6
473  END IF
474  GO TO 10
475  END IF
476 *
477 * Apply permutation to VSL and VSR
478 *
479  IF( ilvsl ) THEN
480  CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
481  $ rwork( iright ), n, vsl, ldvsl, iinfo )
482  IF( iinfo.NE.0 ) THEN
483  info = n + 7
484  GO TO 10
485  END IF
486  END IF
487  IF( ilvsr ) THEN
488  CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
489  $ rwork( iright ), n, vsr, ldvsr, iinfo )
490  IF( iinfo.NE.0 ) THEN
491  info = n + 8
492  GO TO 10
493  END IF
494  END IF
495 *
496 * Undo scaling
497 *
498  IF( ilascl ) THEN
499  CALL zlascl( 'U', -1, -1, anrmto, anrm, n, n, a, lda, iinfo )
500  IF( iinfo.NE.0 ) THEN
501  info = n + 9
502  RETURN
503  END IF
504  CALL zlascl( 'G', -1, -1, anrmto, anrm, n, 1, alpha, n, iinfo )
505  IF( iinfo.NE.0 ) THEN
506  info = n + 9
507  RETURN
508  END IF
509  END IF
510 *
511  IF( ilbscl ) THEN
512  CALL zlascl( 'U', -1, -1, bnrmto, bnrm, n, n, b, ldb, iinfo )
513  IF( iinfo.NE.0 ) THEN
514  info = n + 9
515  RETURN
516  END IF
517  CALL zlascl( 'G', -1, -1, bnrmto, bnrm, n, 1, beta, n, iinfo )
518  IF( iinfo.NE.0 ) THEN
519  info = n + 9
520  RETURN
521  END IF
522  END IF
523 *
524  10 CONTINUE
525  work( 1 ) = lwkopt
526 *
527  RETURN
528 *
529 * End of ZGEGS
530 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

subroutine zgegv ( character  JOBVL,
character  JOBVR,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 This routine is deprecated and has been replaced by routine ZGGEV.

 ZGEGV computes the eigenvalues and, optionally, the left and/or right
 eigenvectors of a complex matrix pair (A,B).
 Given two square matrices A and B,
 the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
 eigenvalues lambda and corresponding (non-zero) eigenvectors x such
 that
    A*x = lambda*B*x.

 An alternate form is to find the eigenvalues mu and corresponding
 eigenvectors y such that
    mu*A*y = B*y.

 These two forms are equivalent with mu = 1/lambda and x = y if
 neither lambda nor mu is zero.  In order to deal with the case that
 lambda or mu is zero or small, two values alpha and beta are returned
 for each eigenvalue, such that lambda = alpha/beta and
 mu = beta/alpha.

 The vectors x and y in the above equations are right eigenvectors of
 the matrix pair (A,B).  Vectors u and v satisfying
    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
 are left eigenvectors of (A,B).

 Note: this routine performs "full balancing" on A and B
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors (returned
                  in VL).
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors (returned
                  in VR).
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the matrix A.
          If JOBVL = 'V' or JOBVR = 'V', then on exit A
          contains the Schur form of A from the generalized Schur
          factorization of the pair (A,B) after balancing.  If no
          eigenvectors were computed, then only the diagonal elements
          of the Schur form will be correct.  See ZGGHRD and ZHGEQZ
          for details.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the matrix B.
          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
          upper triangular matrix obtained from B in the generalized
          Schur factorization of the pair (A,B) after balancing.
          If no eigenvectors were computed, then only the diagonal
          elements of B will be correct.  See ZGGHRD and ZHGEQZ for
          details.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
          The complex scalars alpha that define the eigenvalues of
          GNEP.
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          The complex scalars beta that define the eigenvalues of GNEP.
          
          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
          represent the j-th eigenvalue of the matrix pair (A,B), in
          one of the forms lambda = alpha/beta or mu = beta/alpha.
          Since either lambda or mu may overflow, they should not,
          in general, be computed.
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored
          in the columns of VL, in the same order as their eigenvalues.
          Each eigenvector is scaled so that its largest component has
          abs(real part) + abs(imag. part) = 1, except for eigenvectors
          corresponding to an eigenvalue with alpha = beta = 0, which
          are set to zero.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors x(j) are stored
          in the columns of VR, in the same order as their eigenvalues.
          Each eigenvector is scaled so that its largest component has
          abs(real part) + abs(imag. part) = 1, except for eigenvectors
          corresponding to an eigenvalue with alpha = beta = 0, which
          are set to zero.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[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).
          For good performance, LWORK must generally be larger.
          To compute the optimal value of LWORK, call ILAENV to get
          blocksizes (for ZGEQRF, ZUNMQR, and ZUNGQR.)  Then compute:
          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and ZUNGQR;
          The optimal LWORK is  MAX( 2*N, N*(NB+1) ).

          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (8*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          =1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHA(j) and BETA(j) should be
                correct for j=INFO+1,...,N.
          > N:  errors that usually indicate LAPACK problems:
                =N+1: error return from ZGGBAL
                =N+2: error return from ZGEQRF
                =N+3: error return from ZUNMQR
                =N+4: error return from ZUNGQR
                =N+5: error return from ZGGHRD
                =N+6: error return from ZHGEQZ (other than failed
                                               iteration)
                =N+7: error return from ZTGEVC
                =N+8: error return from ZGGBAK (computing VL)
                =N+9: error return from ZGGBAK (computing VR)
                =N+10: error return from ZLASCL (various calls)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  Balancing
  ---------

  This driver calls ZGGBAL to both permute and scale rows and columns
  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
  and PL*B*R will be upper triangular except for the diagonal blocks
  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
  possible.  The diagonal scaling matrices DL and DR are chosen so
  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
  one (except for the elements that start out zero.)

  After the eigenvalues and eigenvectors of the balanced matrices
  have been computed, ZGGBAK transforms the eigenvectors back to what
  they would have been (in perfect arithmetic) if they had not been
  balanced.

  Contents of A and B on Exit
  -------- -- - --- - -- ----

  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
  both), then on exit the arrays A and B will contain the complex Schur
  form[*] of the "balanced" versions of A and B.  If no eigenvectors
  are computed, then only the diagonal blocks will be correct.

  [*] In other words, upper triangular form.

Definition at line 284 of file zgegv.f.

284 *
285 * -- LAPACK driver routine (version 3.4.0) --
286 * -- LAPACK is a software package provided by Univ. of Tennessee, --
287 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288 * November 2011
289 *
290 * .. Scalar Arguments ..
291  CHARACTER jobvl, jobvr
292  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
293 * ..
294 * .. Array Arguments ..
295  DOUBLE PRECISION rwork( * )
296  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
297  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
298  $ work( * )
299 * ..
300 *
301 * =====================================================================
302 *
303 * .. Parameters ..
304  DOUBLE PRECISION zero, one
305  parameter( zero = 0.0d0, one = 1.0d0 )
306  COMPLEX*16 czero, cone
307  parameter( czero = ( 0.0d0, 0.0d0 ),
308  $ cone = ( 1.0d0, 0.0d0 ) )
309 * ..
310 * .. Local Scalars ..
311  LOGICAL ilimit, ilv, ilvl, ilvr, lquery
312  CHARACTER chtemp
313  INTEGER icols, ihi, iinfo, ijobvl, ijobvr, ileft, ilo,
314  $ in, iright, irows, irwork, itau, iwork, jc, jr,
315  $ lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3
316  DOUBLE PRECISION absai, absar, absb, anrm, anrm1, anrm2, bnrm,
317  $ bnrm1, bnrm2, eps, safmax, safmin, salfai,
318  $ salfar, sbeta, scale, temp
319  COMPLEX*16 x
320 * ..
321 * .. Local Arrays ..
322  LOGICAL ldumma( 1 )
323 * ..
324 * .. External Subroutines ..
325  EXTERNAL xerbla, zgeqrf, zggbak, zggbal, zgghrd, zhgeqz,
327 * ..
328 * .. External Functions ..
329  LOGICAL lsame
330  INTEGER ilaenv
331  DOUBLE PRECISION dlamch, zlange
332  EXTERNAL lsame, ilaenv, dlamch, zlange
333 * ..
334 * .. Intrinsic Functions ..
335  INTRINSIC abs, dble, dcmplx, dimag, int, max
336 * ..
337 * .. Statement Functions ..
338  DOUBLE PRECISION abs1
339 * ..
340 * .. Statement Function definitions ..
341  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
342 * ..
343 * .. Executable Statements ..
344 *
345 * Decode the input arguments
346 *
347  IF( lsame( jobvl, 'N' ) ) THEN
348  ijobvl = 1
349  ilvl = .false.
350  ELSE IF( lsame( jobvl, 'V' ) ) THEN
351  ijobvl = 2
352  ilvl = .true.
353  ELSE
354  ijobvl = -1
355  ilvl = .false.
356  END IF
357 *
358  IF( lsame( jobvr, 'N' ) ) THEN
359  ijobvr = 1
360  ilvr = .false.
361  ELSE IF( lsame( jobvr, 'V' ) ) THEN
362  ijobvr = 2
363  ilvr = .true.
364  ELSE
365  ijobvr = -1
366  ilvr = .false.
367  END IF
368  ilv = ilvl .OR. ilvr
369 *
370 * Test the input arguments
371 *
372  lwkmin = max( 2*n, 1 )
373  lwkopt = lwkmin
374  work( 1 ) = lwkopt
375  lquery = ( lwork.EQ.-1 )
376  info = 0
377  IF( ijobvl.LE.0 ) THEN
378  info = -1
379  ELSE IF( ijobvr.LE.0 ) THEN
380  info = -2
381  ELSE IF( n.LT.0 ) THEN
382  info = -3
383  ELSE IF( lda.LT.max( 1, n ) ) THEN
384  info = -5
385  ELSE IF( ldb.LT.max( 1, n ) ) THEN
386  info = -7
387  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
388  info = -11
389  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
390  info = -13
391  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
392  info = -15
393  END IF
394 *
395  IF( info.EQ.0 ) THEN
396  nb1 = ilaenv( 1, 'ZGEQRF', ' ', n, n, -1, -1 )
397  nb2 = ilaenv( 1, 'ZUNMQR', ' ', n, n, n, -1 )
398  nb3 = ilaenv( 1, 'ZUNGQR', ' ', n, n, n, -1 )
399  nb = max( nb1, nb2, nb3 )
400  lopt = max( 2*n, n*( nb+1 ) )
401  work( 1 ) = lopt
402  END IF
403 *
404  IF( info.NE.0 ) THEN
405  CALL xerbla( 'ZGEGV ', -info )
406  RETURN
407  ELSE IF( lquery ) THEN
408  RETURN
409  END IF
410 *
411 * Quick return if possible
412 *
413  IF( n.EQ.0 )
414  $ RETURN
415 *
416 * Get machine constants
417 *
418  eps = dlamch( 'E' )*dlamch( 'B' )
419  safmin = dlamch( 'S' )
420  safmin = safmin + safmin
421  safmax = one / safmin
422 *
423 * Scale A
424 *
425  anrm = zlange( 'M', n, n, a, lda, rwork )
426  anrm1 = anrm
427  anrm2 = one
428  IF( anrm.LT.one ) THEN
429  IF( safmax*anrm.LT.one ) THEN
430  anrm1 = safmin
431  anrm2 = safmax*anrm
432  END IF
433  END IF
434 *
435  IF( anrm.GT.zero ) THEN
436  CALL zlascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
437  IF( iinfo.NE.0 ) THEN
438  info = n + 10
439  RETURN
440  END IF
441  END IF
442 *
443 * Scale B
444 *
445  bnrm = zlange( 'M', n, n, b, ldb, rwork )
446  bnrm1 = bnrm
447  bnrm2 = one
448  IF( bnrm.LT.one ) THEN
449  IF( safmax*bnrm.LT.one ) THEN
450  bnrm1 = safmin
451  bnrm2 = safmax*bnrm
452  END IF
453  END IF
454 *
455  IF( bnrm.GT.zero ) THEN
456  CALL zlascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
457  IF( iinfo.NE.0 ) THEN
458  info = n + 10
459  RETURN
460  END IF
461  END IF
462 *
463 * Permute the matrix to make it more nearly triangular
464 * Also "balance" the matrix.
465 *
466  ileft = 1
467  iright = n + 1
468  irwork = iright + n
469  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
470  $ rwork( iright ), rwork( irwork ), iinfo )
471  IF( iinfo.NE.0 ) THEN
472  info = n + 1
473  GO TO 80
474  END IF
475 *
476 * Reduce B to triangular form, and initialize VL and/or VR
477 *
478  irows = ihi + 1 - ilo
479  IF( ilv ) THEN
480  icols = n + 1 - ilo
481  ELSE
482  icols = irows
483  END IF
484  itau = 1
485  iwork = itau + irows
486  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
487  $ work( iwork ), lwork+1-iwork, iinfo )
488  IF( iinfo.GE.0 )
489  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
490  IF( iinfo.NE.0 ) THEN
491  info = n + 2
492  GO TO 80
493  END IF
494 *
495  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
496  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
497  $ lwork+1-iwork, iinfo )
498  IF( iinfo.GE.0 )
499  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
500  IF( iinfo.NE.0 ) THEN
501  info = n + 3
502  GO TO 80
503  END IF
504 *
505  IF( ilvl ) THEN
506  CALL zlaset( 'Full', n, n, czero, cone, vl, ldvl )
507  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
508  $ vl( ilo+1, ilo ), ldvl )
509  CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
510  $ work( itau ), work( iwork ), lwork+1-iwork,
511  $ iinfo )
512  IF( iinfo.GE.0 )
513  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
514  IF( iinfo.NE.0 ) THEN
515  info = n + 4
516  GO TO 80
517  END IF
518  END IF
519 *
520  IF( ilvr )
521  $ CALL zlaset( 'Full', n, n, czero, cone, vr, ldvr )
522 *
523 * Reduce to generalized Hessenberg form
524 *
525  IF( ilv ) THEN
526 *
527 * Eigenvectors requested -- work on whole matrix.
528 *
529  CALL zgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
530  $ ldvl, vr, ldvr, iinfo )
531  ELSE
532  CALL zgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
533  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
534  END IF
535  IF( iinfo.NE.0 ) THEN
536  info = n + 5
537  GO TO 80
538  END IF
539 *
540 * Perform QZ algorithm
541 *
542  iwork = itau
543  IF( ilv ) THEN
544  chtemp = 'S'
545  ELSE
546  chtemp = 'E'
547  END IF
548  CALL zhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
549  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwork ),
550  $ lwork+1-iwork, rwork( irwork ), iinfo )
551  IF( iinfo.GE.0 )
552  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
553  IF( iinfo.NE.0 ) THEN
554  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
555  info = iinfo
556  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
557  info = iinfo - n
558  ELSE
559  info = n + 6
560  END IF
561  GO TO 80
562  END IF
563 *
564  IF( ilv ) THEN
565 *
566 * Compute Eigenvectors
567 *
568  IF( ilvl ) THEN
569  IF( ilvr ) THEN
570  chtemp = 'B'
571  ELSE
572  chtemp = 'L'
573  END IF
574  ELSE
575  chtemp = 'R'
576  END IF
577 *
578  CALL ztgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
579  $ vr, ldvr, n, in, work( iwork ), rwork( irwork ),
580  $ iinfo )
581  IF( iinfo.NE.0 ) THEN
582  info = n + 7
583  GO TO 80
584  END IF
585 *
586 * Undo balancing on VL and VR, rescale
587 *
588  IF( ilvl ) THEN
589  CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
590  $ rwork( iright ), n, vl, ldvl, iinfo )
591  IF( iinfo.NE.0 ) THEN
592  info = n + 8
593  GO TO 80
594  END IF
595  DO 30 jc = 1, n
596  temp = zero
597  DO 10 jr = 1, n
598  temp = max( temp, abs1( vl( jr, jc ) ) )
599  10 CONTINUE
600  IF( temp.LT.safmin )
601  $ GO TO 30
602  temp = one / temp
603  DO 20 jr = 1, n
604  vl( jr, jc ) = vl( jr, jc )*temp
605  20 CONTINUE
606  30 CONTINUE
607  END IF
608  IF( ilvr ) THEN
609  CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
610  $ rwork( iright ), n, vr, ldvr, iinfo )
611  IF( iinfo.NE.0 ) THEN
612  info = n + 9
613  GO TO 80
614  END IF
615  DO 60 jc = 1, n
616  temp = zero
617  DO 40 jr = 1, n
618  temp = max( temp, abs1( vr( jr, jc ) ) )
619  40 CONTINUE
620  IF( temp.LT.safmin )
621  $ GO TO 60
622  temp = one / temp
623  DO 50 jr = 1, n
624  vr( jr, jc ) = vr( jr, jc )*temp
625  50 CONTINUE
626  60 CONTINUE
627  END IF
628 *
629 * End of eigenvector calculation
630 *
631  END IF
632 *
633 * Undo scaling in alpha, beta
634 *
635 * Note: this does not give the alpha and beta for the unscaled
636 * problem.
637 *
638 * Un-scaling is limited to avoid underflow in alpha and beta
639 * if they are significant.
640 *
641  DO 70 jc = 1, n
642  absar = abs( dble( alpha( jc ) ) )
643  absai = abs( dimag( alpha( jc ) ) )
644  absb = abs( dble( beta( jc ) ) )
645  salfar = anrm*dble( alpha( jc ) )
646  salfai = anrm*dimag( alpha( jc ) )
647  sbeta = bnrm*dble( beta( jc ) )
648  ilimit = .false.
649  scale = one
650 *
651 * Check for significant underflow in imaginary part of ALPHA
652 *
653  IF( abs( salfai ).LT.safmin .AND. absai.GE.
654  $ max( safmin, eps*absar, eps*absb ) ) THEN
655  ilimit = .true.
656  scale = ( safmin / anrm1 ) / max( safmin, anrm2*absai )
657  END IF
658 *
659 * Check for significant underflow in real part of ALPHA
660 *
661  IF( abs( salfar ).LT.safmin .AND. absar.GE.
662  $ max( safmin, eps*absai, eps*absb ) ) THEN
663  ilimit = .true.
664  scale = max( scale, ( safmin / anrm1 ) /
665  $ max( safmin, anrm2*absar ) )
666  END IF
667 *
668 * Check for significant underflow in BETA
669 *
670  IF( abs( sbeta ).LT.safmin .AND. absb.GE.
671  $ max( safmin, eps*absar, eps*absai ) ) THEN
672  ilimit = .true.
673  scale = max( scale, ( safmin / bnrm1 ) /
674  $ max( safmin, bnrm2*absb ) )
675  END IF
676 *
677 * Check for possible overflow when limiting scaling
678 *
679  IF( ilimit ) THEN
680  temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
681  $ abs( sbeta ) )
682  IF( temp.GT.one )
683  $ scale = scale / temp
684  IF( scale.LT.one )
685  $ ilimit = .false.
686  END IF
687 *
688 * Recompute un-scaled ALPHA, BETA if necessary.
689 *
690  IF( ilimit ) THEN
691  salfar = ( scale*dble( alpha( jc ) ) )*anrm
692  salfai = ( scale*dimag( alpha( jc ) ) )*anrm
693  sbeta = ( scale*beta( jc ) )*bnrm
694  END IF
695  alpha( jc ) = dcmplx( salfar, salfai )
696  beta( jc ) = sbeta
697  70 CONTINUE
698 *
699  80 CONTINUE
700  work( 1 ) = lwkopt
701 *
702  RETURN
703 *
704 * End of ZGEGV
705 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
Definition: ztgevc.f:221
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

subroutine zgges ( character  JOBVSL,
character  JOBVSR,
character  SORT,
external  SELCTG,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  SDIM,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvsl, * )  VSL,
integer  LDVSL,
complex*16, dimension( ldvsr, * )  VSR,
integer  LDVSR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

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

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

Purpose:
 ZGGES computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, the generalized complex Schur
 form (S, T), and optionally left and/or right Schur vectors (VSL
 and VSR). This gives the generalized Schur factorization

         (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )

 where (VSR)**H is the conjugate-transpose of VSR.

 Optionally, it also orders the eigenvalues so that a selected cluster
 of eigenvalues appears in the leading diagonal blocks of the upper
 triangular matrix S and the upper triangular matrix T. The leading
 columns of VSL and VSR then form an unitary basis for the
 corresponding left and right eigenspaces (deflating subspaces).

 (If only the generalized eigenvalues are needed, use the driver
 ZGGEV instead, which is faster.)

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
 or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
 usually represented as the pair (alpha,beta), as there is a
 reasonable interpretation for beta=0, and even for both being zero.

 A pair of matrices (S,T) is in generalized complex Schur form if S
 and T are upper triangular and, in addition, the diagonal elements
 of T are non-negative real numbers.
Parameters
[in]JOBVSL
          JOBVSL is CHARACTER*1
          = 'N':  do not compute the left Schur vectors;
          = 'V':  compute the left Schur vectors.
[in]JOBVSR
          JOBVSR is CHARACTER*1
          = 'N':  do not compute the right Schur vectors;
          = 'V':  compute the right Schur vectors.
[in]SORT
          SORT is CHARACTER*1
          Specifies whether or not to order the eigenvalues on the
          diagonal of the generalized Schur form.
          = 'N':  Eigenvalues are not ordered;
          = 'S':  Eigenvalues are ordered (see SELCTG).
[in]SELCTG
          SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments
          SELCTG must be declared EXTERNAL in the calling subroutine.
          If SORT = 'N', SELCTG is not referenced.
          If SORT = 'S', SELCTG is used to select eigenvalues to sort
          to the top left of the Schur form.
          An eigenvalue ALPHA(j)/BETA(j) is selected if
          SELCTG(ALPHA(j),BETA(j)) is true.

          Note that a selected complex eigenvalue may no longer satisfy
          SELCTG(ALPHA(j),BETA(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 matrices A, B, VSL, and VSR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the first of the pair of matrices.
          On exit, A has been overwritten by its generalized Schur
          form S.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the second of the pair of matrices.
          On exit, B has been overwritten by its generalized Schur
          form T.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]SDIM
          SDIM is INTEGER
          If SORT = 'N', SDIM = 0.
          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
          for which SELCTG is true.
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit,  ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.  ALPHA(j), j=1,...,N  and  BETA(j),
          j=1,...,N  are the diagonals of the complex Schur form (A,B)
          output by ZGGES. The  BETA(j) will be non-negative real.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VSL
          VSL is COMPLEX*16 array, dimension (LDVSL,N)
          If JOBVSL = 'V', VSL will contain the left Schur vectors.
          Not referenced if JOBVSL = 'N'.
[in]LDVSL
          LDVSL is INTEGER
          The leading dimension of the matrix VSL. LDVSL >= 1, and
          if JOBVSL = 'V', LDVSL >= N.
[out]VSR
          VSR is COMPLEX*16 array, dimension (LDVSR,N)
          If JOBVSR = 'V', VSR will contain the right Schur vectors.
          Not referenced if JOBVSR = 'N'.
[in]LDVSR
          LDVSR is INTEGER
          The leading dimension of the matrix VSR. LDVSR >= 1, and
          if JOBVSR = 'V', LDVSR >= N.
[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).
          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (8*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.
          =1,...,N:
                The QZ iteration failed.  (A,B) are not in Schur
                form, but ALPHA(j) and BETA(j) should be correct for
                j=INFO+1,...,N.
          > N:  =N+1: other than QZ iteration failed in ZHGEQZ
                =N+2: after reordering, roundoff changed values of
                      some complex eigenvalues so that leading
                      eigenvalues in the Generalized Schur form no
                      longer satisfy SELCTG=.TRUE.  This could also
                      be caused due to scaling.
                =N+3: reordering failed in ZTGSEN.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 272 of file zgges.f.

272 *
273 * -- LAPACK driver routine (version 3.6.0) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * November 2015
277 *
278 * .. Scalar Arguments ..
279  CHARACTER jobvsl, jobvsr, sort
280  INTEGER info, lda, ldb, ldvsl, ldvsr, lwork, n, sdim
281 * ..
282 * .. Array Arguments ..
283  LOGICAL bwork( * )
284  DOUBLE PRECISION rwork( * )
285  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
286  $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
287  $ work( * )
288 * ..
289 * .. Function Arguments ..
290  LOGICAL selctg
291  EXTERNAL selctg
292 * ..
293 *
294 * =====================================================================
295 *
296 * .. Parameters ..
297  DOUBLE PRECISION zero, one
298  parameter( zero = 0.0d0, one = 1.0d0 )
299  COMPLEX*16 czero, cone
300  parameter( czero = ( 0.0d0, 0.0d0 ),
301  $ cone = ( 1.0d0, 0.0d0 ) )
302 * ..
303 * .. Local Scalars ..
304  LOGICAL cursl, ilascl, ilbscl, ilvsl, ilvsr, lastsl,
305  $ lquery, wantst
306  INTEGER i, icols, ierr, ihi, ijobvl, ijobvr, ileft,
307  $ ilo, iright, irows, irwrk, itau, iwrk, lwkmin,
308  $ lwkopt
309  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps, pvsl,
310  $ pvsr, smlnum
311 * ..
312 * .. Local Arrays ..
313  INTEGER idum( 1 )
314  DOUBLE PRECISION dif( 2 )
315 * ..
316 * .. External Subroutines ..
317  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghrd,
319  $ zunmqr
320 * ..
321 * .. External Functions ..
322  LOGICAL lsame
323  INTEGER ilaenv
324  DOUBLE PRECISION dlamch, zlange
325  EXTERNAL lsame, ilaenv, dlamch, zlange
326 * ..
327 * .. Intrinsic Functions ..
328  INTRINSIC max, sqrt
329 * ..
330 * .. Executable Statements ..
331 *
332 * Decode the input arguments
333 *
334  IF( lsame( jobvsl, 'N' ) ) THEN
335  ijobvl = 1
336  ilvsl = .false.
337  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
338  ijobvl = 2
339  ilvsl = .true.
340  ELSE
341  ijobvl = -1
342  ilvsl = .false.
343  END IF
344 *
345  IF( lsame( jobvsr, 'N' ) ) THEN
346  ijobvr = 1
347  ilvsr = .false.
348  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
349  ijobvr = 2
350  ilvsr = .true.
351  ELSE
352  ijobvr = -1
353  ilvsr = .false.
354  END IF
355 *
356  wantst = lsame( sort, 'S' )
357 *
358 * Test the input arguments
359 *
360  info = 0
361  lquery = ( lwork.EQ.-1 )
362  IF( ijobvl.LE.0 ) THEN
363  info = -1
364  ELSE IF( ijobvr.LE.0 ) THEN
365  info = -2
366  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
367  info = -3
368  ELSE IF( n.LT.0 ) THEN
369  info = -5
370  ELSE IF( lda.LT.max( 1, n ) ) THEN
371  info = -7
372  ELSE IF( ldb.LT.max( 1, n ) ) THEN
373  info = -9
374  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
375  info = -14
376  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
377  info = -16
378  END IF
379 *
380 * Compute workspace
381 * (Note: Comments in the code beginning "Workspace:" describe the
382 * minimal amount of workspace needed at that point in the code,
383 * as well as the preferred amount for good performance.
384 * NB refers to the optimal block size for the immediately
385 * following subroutine, as returned by ILAENV.)
386 *
387  IF( info.EQ.0 ) THEN
388  lwkmin = max( 1, 2*n )
389  lwkopt = max( 1, n + n*ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
390  lwkopt = max( lwkopt, n +
391  $ n*ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, -1 ) )
392  IF( ilvsl ) THEN
393  lwkopt = max( lwkopt, n +
394  $ n*ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, -1 ) )
395  END IF
396  work( 1 ) = lwkopt
397 *
398  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
399  $ info = -18
400  END IF
401 *
402  IF( info.NE.0 ) THEN
403  CALL xerbla( 'ZGGES ', -info )
404  RETURN
405  ELSE IF( lquery ) THEN
406  RETURN
407  END IF
408 *
409 * Quick return if possible
410 *
411  IF( n.EQ.0 ) THEN
412  sdim = 0
413  RETURN
414  END IF
415 *
416 * Get machine constants
417 *
418  eps = dlamch( 'P' )
419  smlnum = dlamch( 'S' )
420  bignum = one / smlnum
421  CALL dlabad( smlnum, bignum )
422  smlnum = sqrt( smlnum ) / eps
423  bignum = one / smlnum
424 *
425 * Scale A if max element outside range [SMLNUM,BIGNUM]
426 *
427  anrm = zlange( 'M', n, n, a, lda, rwork )
428  ilascl = .false.
429  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
430  anrmto = smlnum
431  ilascl = .true.
432  ELSE IF( anrm.GT.bignum ) THEN
433  anrmto = bignum
434  ilascl = .true.
435  END IF
436 *
437  IF( ilascl )
438  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
439 *
440 * Scale B if max element outside range [SMLNUM,BIGNUM]
441 *
442  bnrm = zlange( 'M', n, n, b, ldb, rwork )
443  ilbscl = .false.
444  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
445  bnrmto = smlnum
446  ilbscl = .true.
447  ELSE IF( bnrm.GT.bignum ) THEN
448  bnrmto = bignum
449  ilbscl = .true.
450  END IF
451 *
452  IF( ilbscl )
453  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
454 *
455 * Permute the matrix to make it more nearly triangular
456 * (Real Workspace: need 6*N)
457 *
458  ileft = 1
459  iright = n + 1
460  irwrk = iright + n
461  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
462  $ rwork( iright ), rwork( irwrk ), ierr )
463 *
464 * Reduce B to triangular form (QR decomposition of B)
465 * (Complex Workspace: need N, prefer N*NB)
466 *
467  irows = ihi + 1 - ilo
468  icols = n + 1 - ilo
469  itau = 1
470  iwrk = itau + irows
471  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
472  $ work( iwrk ), lwork+1-iwrk, ierr )
473 *
474 * Apply the orthogonal transformation to matrix A
475 * (Complex Workspace: need N, prefer N*NB)
476 *
477  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
478  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
479  $ lwork+1-iwrk, ierr )
480 *
481 * Initialize VSL
482 * (Complex Workspace: need N, prefer N*NB)
483 *
484  IF( ilvsl ) THEN
485  CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
486  IF( irows.GT.1 ) THEN
487  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
488  $ vsl( ilo+1, ilo ), ldvsl )
489  END IF
490  CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
491  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
492  END IF
493 *
494 * Initialize VSR
495 *
496  IF( ilvsr )
497  $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
498 *
499 * Reduce to generalized Hessenberg form
500 * (Workspace: none needed)
501 *
502  CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
503  $ ldvsl, vsr, ldvsr, ierr )
504 *
505  sdim = 0
506 *
507 * Perform QZ algorithm, computing Schur vectors if desired
508 * (Complex Workspace: need N)
509 * (Real Workspace: need N)
510 *
511  iwrk = itau
512  CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
513  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
514  $ lwork+1-iwrk, rwork( irwrk ), ierr )
515  IF( ierr.NE.0 ) THEN
516  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
517  info = ierr
518  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
519  info = ierr - n
520  ELSE
521  info = n + 1
522  END IF
523  GO TO 30
524  END IF
525 *
526 * Sort eigenvalues ALPHA/BETA if desired
527 * (Workspace: none needed)
528 *
529  IF( wantst ) THEN
530 *
531 * Undo scaling on eigenvalues before selecting
532 *
533  IF( ilascl )
534  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
535  IF( ilbscl )
536  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, 1, beta, n, ierr )
537 *
538 * Select eigenvalues
539 *
540  DO 10 i = 1, n
541  bwork( i ) = selctg( alpha( i ), beta( i ) )
542  10 CONTINUE
543 *
544  CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alpha,
545  $ beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl, pvsr,
546  $ dif, work( iwrk ), lwork-iwrk+1, idum, 1, ierr )
547  IF( ierr.EQ.1 )
548  $ info = n + 3
549 *
550  END IF
551 *
552 * Apply back-permutation to VSL and VSR
553 * (Workspace: none needed)
554 *
555  IF( ilvsl )
556  $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
557  $ rwork( iright ), n, vsl, ldvsl, ierr )
558  IF( ilvsr )
559  $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
560  $ rwork( iright ), n, vsr, ldvsr, ierr )
561 *
562 * Undo scaling
563 *
564  IF( ilascl ) THEN
565  CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
566  CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
567  END IF
568 *
569  IF( ilbscl ) THEN
570  CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
571  CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
572  END IF
573 *
574  IF( wantst ) THEN
575 *
576 * Check if reordering is correct
577 *
578  lastsl = .true.
579  sdim = 0
580  DO 20 i = 1, n
581  cursl = selctg( alpha( i ), beta( i ) )
582  IF( cursl )
583  $ sdim = sdim + 1
584  IF( cursl .AND. .NOT.lastsl )
585  $ info = n + 2
586  lastsl = cursl
587  20 CONTINUE
588 *
589  END IF
590 *
591  30 CONTINUE
592 *
593  work( 1 ) = lwkopt
594 *
595  RETURN
596 *
597 * End of ZGGES
598 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine ztgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
ZTGSEN
Definition: ztgsen.f:435
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgges3 ( character  JOBVSL,
character  JOBVSR,
character  SORT,
external  SELCTG,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  SDIM,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvsl, * )  VSL,
integer  LDVSL,
complex*16, dimension( ldvsr, * )  VSR,
integer  LDVSR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

ZGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)

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

Purpose:
 ZGGES3 computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, the generalized complex Schur
 form (S, T), and optionally left and/or right Schur vectors (VSL
 and VSR). This gives the generalized Schur factorization

         (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )

 where (VSR)**H is the conjugate-transpose of VSR.

 Optionally, it also orders the eigenvalues so that a selected cluster
 of eigenvalues appears in the leading diagonal blocks of the upper
 triangular matrix S and the upper triangular matrix T. The leading
 columns of VSL and VSR then form an unitary basis for the
 corresponding left and right eigenspaces (deflating subspaces).

 (If only the generalized eigenvalues are needed, use the driver
 ZGGEV instead, which is faster.)

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
 or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
 usually represented as the pair (alpha,beta), as there is a
 reasonable interpretation for beta=0, and even for both being zero.

 A pair of matrices (S,T) is in generalized complex Schur form if S
 and T are upper triangular and, in addition, the diagonal elements
 of T are non-negative real numbers.
Parameters
[in]JOBVSL
          JOBVSL is CHARACTER*1
          = 'N':  do not compute the left Schur vectors;
          = 'V':  compute the left Schur vectors.
[in]JOBVSR
          JOBVSR is CHARACTER*1
          = 'N':  do not compute the right Schur vectors;
          = 'V':  compute the right Schur vectors.
[in]SORT
          SORT is CHARACTER*1
          Specifies whether or not to order the eigenvalues on the
          diagonal of the generalized Schur form.
          = 'N':  Eigenvalues are not ordered;
          = 'S':  Eigenvalues are ordered (see SELCTG).
[in]SELCTG
          SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments
          SELCTG must be declared EXTERNAL in the calling subroutine.
          If SORT = 'N', SELCTG is not referenced.
          If SORT = 'S', SELCTG is used to select eigenvalues to sort
          to the top left of the Schur form.
          An eigenvalue ALPHA(j)/BETA(j) is selected if
          SELCTG(ALPHA(j),BETA(j)) is true.

          Note that a selected complex eigenvalue may no longer satisfy
          SELCTG(ALPHA(j),BETA(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 matrices A, B, VSL, and VSR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the first of the pair of matrices.
          On exit, A has been overwritten by its generalized Schur
          form S.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the second of the pair of matrices.
          On exit, B has been overwritten by its generalized Schur
          form T.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]SDIM
          SDIM is INTEGER
          If SORT = 'N', SDIM = 0.
          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
          for which SELCTG is true.
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit,  ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.  ALPHA(j), j=1,...,N  and  BETA(j),
          j=1,...,N  are the diagonals of the complex Schur form (A,B)
          output by ZGGES3. The  BETA(j) will be non-negative real.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VSL
          VSL is COMPLEX*16 array, dimension (LDVSL,N)
          If JOBVSL = 'V', VSL will contain the left Schur vectors.
          Not referenced if JOBVSL = 'N'.
[in]LDVSL
          LDVSL is INTEGER
          The leading dimension of the matrix VSL. LDVSL >= 1, and
          if JOBVSL = 'V', LDVSL >= N.
[out]VSR
          VSR is COMPLEX*16 array, dimension (LDVSR,N)
          If JOBVSR = 'V', VSR will contain the right Schur vectors.
          Not referenced if JOBVSR = 'N'.
[in]LDVSR
          LDVSR is INTEGER
          The leading dimension of the matrix VSR. LDVSR >= 1, and
          if JOBVSR = 'V', LDVSR >= N.
[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.

          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (8*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.
          =1,...,N:
                The QZ iteration failed.  (A,B) are not in Schur
                form, but ALPHA(j) and BETA(j) should be correct for
                j=INFO+1,...,N.
          > N:  =N+1: other than QZ iteration failed in ZHGEQZ
                =N+2: after reordering, roundoff changed values of
                      some complex eigenvalues so that leading
                      eigenvalues in the Generalized Schur form no
                      longer satisfy SELCTG=.TRUE.  This could also
                      be caused due to scaling.
                =N+3: reordering failed in ZTGSEN.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
January 2015

Definition at line 271 of file zgges3.f.

271 *
272 * -- LAPACK driver routine (version 3.6.0) --
273 * -- LAPACK is a software package provided by Univ. of Tennessee, --
274 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
275 * January 2015
276 *
277 * .. Scalar Arguments ..
278  CHARACTER jobvsl, jobvsr, sort
279  INTEGER info, lda, ldb, ldvsl, ldvsr, lwork, n, sdim
280 * ..
281 * .. Array Arguments ..
282  LOGICAL bwork( * )
283  DOUBLE PRECISION rwork( * )
284  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
285  $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
286  $ work( * )
287 * ..
288 * .. Function Arguments ..
289  LOGICAL selctg
290  EXTERNAL selctg
291 * ..
292 *
293 * =====================================================================
294 *
295 * .. Parameters ..
296  DOUBLE PRECISION zero, one
297  parameter( zero = 0.0d0, one = 1.0d0 )
298  COMPLEX*16 czero, cone
299  parameter( czero = ( 0.0d0, 0.0d0 ),
300  $ cone = ( 1.0d0, 0.0d0 ) )
301 * ..
302 * .. Local Scalars ..
303  LOGICAL cursl, ilascl, ilbscl, ilvsl, ilvsr, lastsl,
304  $ lquery, wantst
305  INTEGER i, icols, ierr, ihi, ijobvl, ijobvr, ileft,
306  $ ilo, iright, irows, irwrk, itau, iwrk, lwkopt
307  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps, pvsl,
308  $ pvsr, smlnum
309 * ..
310 * .. Local Arrays ..
311  INTEGER idum( 1 )
312  DOUBLE PRECISION dif( 2 )
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghd3,
317  $ zunmqr
318 * ..
319 * .. External Functions ..
320  LOGICAL lsame
321  DOUBLE PRECISION dlamch, zlange
322  EXTERNAL lsame, dlamch, zlange
323 * ..
324 * .. Intrinsic Functions ..
325  INTRINSIC max, sqrt
326 * ..
327 * .. Executable Statements ..
328 *
329 * Decode the input arguments
330 *
331  IF( lsame( jobvsl, 'N' ) ) THEN
332  ijobvl = 1
333  ilvsl = .false.
334  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
335  ijobvl = 2
336  ilvsl = .true.
337  ELSE
338  ijobvl = -1
339  ilvsl = .false.
340  END IF
341 *
342  IF( lsame( jobvsr, 'N' ) ) THEN
343  ijobvr = 1
344  ilvsr = .false.
345  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
346  ijobvr = 2
347  ilvsr = .true.
348  ELSE
349  ijobvr = -1
350  ilvsr = .false.
351  END IF
352 *
353  wantst = lsame( sort, 'S' )
354 *
355 * Test the input arguments
356 *
357  info = 0
358  lquery = ( lwork.EQ.-1 )
359  IF( ijobvl.LE.0 ) THEN
360  info = -1
361  ELSE IF( ijobvr.LE.0 ) THEN
362  info = -2
363  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
364  info = -3
365  ELSE IF( n.LT.0 ) THEN
366  info = -5
367  ELSE IF( lda.LT.max( 1, n ) ) THEN
368  info = -7
369  ELSE IF( ldb.LT.max( 1, n ) ) THEN
370  info = -9
371  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
372  info = -14
373  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
374  info = -16
375  ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
376  info = -18
377  END IF
378 *
379 * Compute workspace
380 *
381  IF( info.EQ.0 ) THEN
382  CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
383  lwkopt = max( 1, n + int( work( 1 ) ) )
384  CALL zunmqr( 'L', 'C', n, n, n, b, ldb, work, a, lda, work,
385  $ -1, ierr )
386  lwkopt = max( lwkopt, n + int( work( 1 ) ) )
387  IF( ilvsl ) THEN
388  CALL zungqr( n, n, n, vsl, ldvsl, work, work, -1, ierr )
389  lwkopt = max( lwkopt, n + int( work( 1 ) ) )
390  END IF
391  CALL zgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
392  $ ldvsl, vsr, ldvsr, work, -1, ierr )
393  lwkopt = max( lwkopt, n + int( work( 1 ) ) )
394  CALL zhgeqz( 'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
395  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work,
396  $ -1, rwork, ierr )
397  lwkopt = max( lwkopt, int( work( 1 ) ) )
398  IF( wantst ) THEN
399  CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
400  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, sdim,
401  $ pvsl, pvsr, dif, work, -1, idum, 1, ierr )
402  lwkopt = max( lwkopt, int( work( 1 ) ) )
403  END IF
404  work( 1 ) = dcmplx( lwkopt )
405  END IF
406 *
407  IF( info.NE.0 ) THEN
408  CALL xerbla( 'ZGGES3 ', -info )
409  RETURN
410  ELSE IF( lquery ) THEN
411  RETURN
412  END IF
413 *
414 * Quick return if possible
415 *
416  IF( n.EQ.0 ) THEN
417  sdim = 0
418  RETURN
419  END IF
420 *
421 * Get machine constants
422 *
423  eps = dlamch( 'P' )
424  smlnum = dlamch( 'S' )
425  bignum = one / smlnum
426  CALL dlabad( smlnum, bignum )
427  smlnum = sqrt( smlnum ) / eps
428  bignum = one / smlnum
429 *
430 * Scale A if max element outside range [SMLNUM,BIGNUM]
431 *
432  anrm = zlange( 'M', n, n, a, lda, rwork )
433  ilascl = .false.
434  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
435  anrmto = smlnum
436  ilascl = .true.
437  ELSE IF( anrm.GT.bignum ) THEN
438  anrmto = bignum
439  ilascl = .true.
440  END IF
441 *
442  IF( ilascl )
443  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
444 *
445 * Scale B if max element outside range [SMLNUM,BIGNUM]
446 *
447  bnrm = zlange( 'M', n, n, b, ldb, rwork )
448  ilbscl = .false.
449  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
450  bnrmto = smlnum
451  ilbscl = .true.
452  ELSE IF( bnrm.GT.bignum ) THEN
453  bnrmto = bignum
454  ilbscl = .true.
455  END IF
456 *
457  IF( ilbscl )
458  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
459 *
460 * Permute the matrix to make it more nearly triangular
461 *
462  ileft = 1
463  iright = n + 1
464  irwrk = iright + n
465  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
466  $ rwork( iright ), rwork( irwrk ), ierr )
467 *
468 * Reduce B to triangular form (QR decomposition of B)
469 *
470  irows = ihi + 1 - ilo
471  icols = n + 1 - ilo
472  itau = 1
473  iwrk = itau + irows
474  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
475  $ work( iwrk ), lwork+1-iwrk, ierr )
476 *
477 * Apply the orthogonal transformation to matrix A
478 *
479  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
480  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
481  $ lwork+1-iwrk, ierr )
482 *
483 * Initialize VSL
484 *
485  IF( ilvsl ) THEN
486  CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
487  IF( irows.GT.1 ) THEN
488  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
489  $ vsl( ilo+1, ilo ), ldvsl )
490  END IF
491  CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
492  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
493  END IF
494 *
495 * Initialize VSR
496 *
497  IF( ilvsr )
498  $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
499 *
500 * Reduce to generalized Hessenberg form
501 *
502  CALL zgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
503  $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk, ierr )
504 *
505  sdim = 0
506 *
507 * Perform QZ algorithm, computing Schur vectors if desired
508 *
509  iwrk = itau
510  CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
511  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
512  $ lwork+1-iwrk, rwork( irwrk ), ierr )
513  IF( ierr.NE.0 ) THEN
514  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
515  info = ierr
516  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
517  info = ierr - n
518  ELSE
519  info = n + 1
520  END IF
521  GO TO 30
522  END IF
523 *
524 * Sort eigenvalues ALPHA/BETA if desired
525 *
526  IF( wantst ) THEN
527 *
528 * Undo scaling on eigenvalues before selecting
529 *
530  IF( ilascl )
531  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
532  IF( ilbscl )
533  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, 1, beta, n, ierr )
534 *
535 * Select eigenvalues
536 *
537  DO 10 i = 1, n
538  bwork( i ) = selctg( alpha( i ), beta( i ) )
539  10 CONTINUE
540 *
541  CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alpha,
542  $ beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl, pvsr,
543  $ dif, work( iwrk ), lwork-iwrk+1, idum, 1, ierr )
544  IF( ierr.EQ.1 )
545  $ info = n + 3
546 *
547  END IF
548 *
549 * Apply back-permutation to VSL and VSR
550 *
551  IF( ilvsl )
552  $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
553  $ rwork( iright ), n, vsl, ldvsl, ierr )
554  IF( ilvsr )
555  $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
556  $ rwork( iright ), n, vsr, ldvsr, ierr )
557 *
558 * Undo scaling
559 *
560  IF( ilascl ) THEN
561  CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
562  CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
563  END IF
564 *
565  IF( ilbscl ) THEN
566  CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
567  CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
568  END IF
569 *
570  IF( wantst ) THEN
571 *
572 * Check if reordering is correct
573 *
574  lastsl = .true.
575  sdim = 0
576  DO 20 i = 1, n
577  cursl = selctg( alpha( i ), beta( i ) )
578  IF( cursl )
579  $ sdim = sdim + 1
580  IF( cursl .AND. .NOT.lastsl )
581  $ info = n + 2
582  lastsl = cursl
583  20 CONTINUE
584 *
585  END IF
586 *
587  30 CONTINUE
588 *
589  work( 1 ) = dcmplx( lwkopt )
590 *
591  RETURN
592 *
593 * End of ZGGES3
594 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine ztgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
ZTGSEN
Definition: ztgsen.f:435
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
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:141
subroutine zgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
ZGGHD3
Definition: zgghd3.f:229
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zggesx ( character  JOBVSL,
character  JOBVSR,
character  SORT,
external  SELCTG,
character  SENSE,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  SDIM,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvsl, * )  VSL,
integer  LDVSL,
complex*16, dimension( ldvsr, * )  VSR,
integer  LDVSR,
double precision, dimension( 2 )  RCONDE,
double precision, dimension( 2 )  RCONDV,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

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

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

Purpose:
 ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, the complex Schur form (S,T),
 and, optionally, the left and/or right matrices of Schur vectors (VSL
 and VSR).  This gives the generalized Schur factorization

      (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )

 where (VSR)**H is the conjugate-transpose of VSR.

 Optionally, it also orders the eigenvalues so that a selected cluster
 of eigenvalues appears in the leading diagonal blocks of the upper
 triangular matrix S and the upper triangular matrix T; computes
 a reciprocal condition number for the average of the selected
 eigenvalues (RCONDE); and computes a reciprocal condition number for
 the right and left deflating subspaces corresponding to the selected
 eigenvalues (RCONDV). The leading columns of VSL and VSR then form
 an orthonormal basis for the corresponding left and right eigenspaces
 (deflating subspaces).

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
 or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
 usually represented as the pair (alpha,beta), as there is a
 reasonable interpretation for beta=0 or for both being zero.

 A pair of matrices (S,T) is in generalized complex Schur form if T is
 upper triangular with non-negative diagonal and S is upper
 triangular.
Parameters
[in]JOBVSL
          JOBVSL is CHARACTER*1
          = 'N':  do not compute the left Schur vectors;
          = 'V':  compute the left Schur vectors.
[in]JOBVSR
          JOBVSR is CHARACTER*1
          = 'N':  do not compute the right Schur vectors;
          = 'V':  compute the right Schur vectors.
[in]SORT
          SORT is CHARACTER*1
          Specifies whether or not to order the eigenvalues on the
          diagonal of the generalized Schur form.
          = 'N':  Eigenvalues are not ordered;
          = 'S':  Eigenvalues are ordered (see SELCTG).
[in]SELCTG
          SELCTG is procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
          SELCTG must be declared EXTERNAL in the calling subroutine.
          If SORT = 'N', SELCTG is not referenced.
          If SORT = 'S', SELCTG is used to select eigenvalues to sort
          to the top left of the Schur form.
          Note that a selected complex eigenvalue may no longer satisfy
          SELCTG(ALPHA(j),BETA(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+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 deflating subspaces only;
          = 'B' : Computed for both.
          If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
[in]N
          N is INTEGER
          The order of the matrices A, B, VSL, and VSR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the first of the pair of matrices.
          On exit, A has been overwritten by its generalized Schur
          form S.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the second of the pair of matrices.
          On exit, B has been overwritten by its generalized Schur
          form T.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]SDIM
          SDIM is INTEGER
          If SORT = 'N', SDIM = 0.
          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
          for which SELCTG is true.
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.  ALPHA(j) and BETA(j),j=1,...,N  are
          the diagonals of the complex Schur form (S,T).  BETA(j) will
          be non-negative real.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VSL
          VSL is COMPLEX*16 array, dimension (LDVSL,N)
          If JOBVSL = 'V', VSL will contain the left Schur vectors.
          Not referenced if JOBVSL = 'N'.
[in]LDVSL
          LDVSL is INTEGER
          The leading dimension of the matrix VSL. LDVSL >=1, and
          if JOBVSL = 'V', LDVSL >= N.
[out]VSR
          VSR is COMPLEX*16 array, dimension (LDVSR,N)
          If JOBVSR = 'V', VSR will contain the right Schur vectors.
          Not referenced if JOBVSR = 'N'.
[in]LDVSR
          LDVSR is INTEGER
          The leading dimension of the matrix VSR. LDVSR >= 1, and
          if JOBVSR = 'V', LDVSR >= N.
[out]RCONDE
          RCONDE is DOUBLE PRECISION array, dimension ( 2 )
          If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
          reciprocal condition numbers for the average of the selected
          eigenvalues.
          Not referenced if SENSE = 'N' or 'V'.
[out]RCONDV
          RCONDV is DOUBLE PRECISION array, dimension ( 2 )
          If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
          reciprocal condition number for the selected deflating
          subspaces.
          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.
          If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
          LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else
          LWORK >= MAX(1,2*N).  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.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the bound on the optimal size of the WORK
          array and the minimum size of the IWORK array, returns these
          values as the first entries of the WORK and IWORK arrays, and
          no error message related to LWORK or LIWORK is issued by
          XERBLA.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension ( 8*N )
          Real workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.
          If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
          LIWORK >= N+2.

          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the bound on the optimal size of the
          WORK array and the minimum size of the IWORK array, returns
          these values as the first entries of the WORK and IWORK
          arrays, and no error message related to LWORK or LIWORK 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.
          = 1,...,N:
                The QZ iteration failed.  (A,B) are not in Schur
                form, but ALPHA(j) and BETA(j) should be correct for
                j=INFO+1,...,N.
          > N:  =N+1: other than QZ iteration failed in ZHGEQZ
                =N+2: after reordering, roundoff changed values of
                      some complex eigenvalues so that leading
                      eigenvalues in the Generalized Schur form no
                      longer satisfy SELCTG=.TRUE.  This could also
                      be caused due to scaling.
                =N+3: reordering failed in ZTGSEN.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 332 of file zggesx.f.

332 *
333 * -- LAPACK driver routine (version 3.4.0) --
334 * -- LAPACK is a software package provided by Univ. of Tennessee, --
335 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
336 * November 2011
337 *
338 * .. Scalar Arguments ..
339  CHARACTER jobvsl, jobvsr, sense, sort
340  INTEGER info, lda, ldb, ldvsl, ldvsr, liwork, lwork, n,
341  $ sdim
342 * ..
343 * .. Array Arguments ..
344  LOGICAL bwork( * )
345  INTEGER iwork( * )
346  DOUBLE PRECISION rconde( 2 ), rcondv( 2 ), rwork( * )
347  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
348  $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
349  $ work( * )
350 * ..
351 * .. Function Arguments ..
352  LOGICAL selctg
353  EXTERNAL selctg
354 * ..
355 *
356 * =====================================================================
357 *
358 * .. Parameters ..
359  DOUBLE PRECISION zero, one
360  parameter( zero = 0.0d+0, one = 1.0d+0 )
361  COMPLEX*16 czero, cone
362  parameter( czero = ( 0.0d+0, 0.0d+0 ),
363  $ cone = ( 1.0d+0, 0.0d+0 ) )
364 * ..
365 * .. Local Scalars ..
366  LOGICAL cursl, ilascl, ilbscl, ilvsl, ilvsr, lastsl,
367  $ lquery, wantsb, wantse, wantsn, wantst, wantsv
368  INTEGER i, icols, ierr, ihi, ijob, ijobvl, ijobvr,
369  $ ileft, ilo, iright, irows, irwrk, itau, iwrk,
370  $ liwmin, lwrk, maxwrk, minwrk
371  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps, pl,
372  $ pr, smlnum
373 * ..
374 * .. Local Arrays ..
375  DOUBLE PRECISION dif( 2 )
376 * ..
377 * .. External Subroutines ..
378  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghrd,
380  $ zunmqr
381 * ..
382 * .. External Functions ..
383  LOGICAL lsame
384  INTEGER ilaenv
385  DOUBLE PRECISION dlamch, zlange
386  EXTERNAL lsame, ilaenv, dlamch, zlange
387 * ..
388 * .. Intrinsic Functions ..
389  INTRINSIC max, sqrt
390 * ..
391 * .. Executable Statements ..
392 *
393 * Decode the input arguments
394 *
395  IF( lsame( jobvsl, 'N' ) ) THEN
396  ijobvl = 1
397  ilvsl = .false.
398  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
399  ijobvl = 2
400  ilvsl = .true.
401  ELSE
402  ijobvl = -1
403  ilvsl = .false.
404  END IF
405 *
406  IF( lsame( jobvsr, 'N' ) ) THEN
407  ijobvr = 1
408  ilvsr = .false.
409  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
410  ijobvr = 2
411  ilvsr = .true.
412  ELSE
413  ijobvr = -1
414  ilvsr = .false.
415  END IF
416 *
417  wantst = lsame( sort, 'S' )
418  wantsn = lsame( sense, 'N' )
419  wantse = lsame( sense, 'E' )
420  wantsv = lsame( sense, 'V' )
421  wantsb = lsame( sense, 'B' )
422  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
423  IF( wantsn ) THEN
424  ijob = 0
425  ELSE IF( wantse ) THEN
426  ijob = 1
427  ELSE IF( wantsv ) THEN
428  ijob = 2
429  ELSE IF( wantsb ) THEN
430  ijob = 4
431  END IF
432 *
433 * Test the input arguments
434 *
435  info = 0
436  IF( ijobvl.LE.0 ) THEN
437  info = -1
438  ELSE IF( ijobvr.LE.0 ) THEN
439  info = -2
440  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
441  info = -3
442  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
443  $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
444  info = -5
445  ELSE IF( n.LT.0 ) THEN
446  info = -6
447  ELSE IF( lda.LT.max( 1, n ) ) THEN
448  info = -8
449  ELSE IF( ldb.LT.max( 1, n ) ) THEN
450  info = -10
451  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
452  info = -15
453  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
454  info = -17
455  END IF
456 *
457 * Compute workspace
458 * (Note: Comments in the code beginning "Workspace:" describe the
459 * minimal amount of workspace needed at that point in the code,
460 * as well as the preferred amount for good performance.
461 * NB refers to the optimal block size for the immediately
462 * following subroutine, as returned by ILAENV.)
463 *
464  IF( info.EQ.0 ) THEN
465  IF( n.GT.0) THEN
466  minwrk = 2*n
467  maxwrk = n*(1 + ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
468  maxwrk = max( maxwrk, n*( 1 +
469  $ ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, -1 ) ) )
470  IF( ilvsl ) THEN
471  maxwrk = max( maxwrk, n*( 1 +
472  $ ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, -1 ) ) )
473  END IF
474  lwrk = maxwrk
475  IF( ijob.GE.1 )
476  $ lwrk = max( lwrk, n*n/2 )
477  ELSE
478  minwrk = 1
479  maxwrk = 1
480  lwrk = 1
481  END IF
482  work( 1 ) = lwrk
483  IF( wantsn .OR. n.EQ.0 ) THEN
484  liwmin = 1
485  ELSE
486  liwmin = n + 2
487  END IF
488  iwork( 1 ) = liwmin
489 *
490  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
491  info = -21
492  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery) THEN
493  info = -24
494  END IF
495  END IF
496 *
497  IF( info.NE.0 ) THEN
498  CALL xerbla( 'ZGGESX', -info )
499  RETURN
500  ELSE IF (lquery) THEN
501  RETURN
502  END IF
503 *
504 * Quick return if possible
505 *
506  IF( n.EQ.0 ) THEN
507  sdim = 0
508  RETURN
509  END IF
510 *
511 * Get machine constants
512 *
513  eps = dlamch( 'P' )
514  smlnum = dlamch( 'S' )
515  bignum = one / smlnum
516  CALL dlabad( smlnum, bignum )
517  smlnum = sqrt( smlnum ) / eps
518  bignum = one / smlnum
519 *
520 * Scale A if max element outside range [SMLNUM,BIGNUM]
521 *
522  anrm = zlange( 'M', n, n, a, lda, rwork )
523  ilascl = .false.
524  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
525  anrmto = smlnum
526  ilascl = .true.
527  ELSE IF( anrm.GT.bignum ) THEN
528  anrmto = bignum
529  ilascl = .true.
530  END IF
531  IF( ilascl )
532  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
533 *
534 * Scale B if max element outside range [SMLNUM,BIGNUM]
535 *
536  bnrm = zlange( 'M', n, n, b, ldb, rwork )
537  ilbscl = .false.
538  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
539  bnrmto = smlnum
540  ilbscl = .true.
541  ELSE IF( bnrm.GT.bignum ) THEN
542  bnrmto = bignum
543  ilbscl = .true.
544  END IF
545  IF( ilbscl )
546  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
547 *
548 * Permute the matrix to make it more nearly triangular
549 * (Real Workspace: need 6*N)
550 *
551  ileft = 1
552  iright = n + 1
553  irwrk = iright + n
554  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
555  $ rwork( iright ), rwork( irwrk ), ierr )
556 *
557 * Reduce B to triangular form (QR decomposition of B)
558 * (Complex Workspace: need N, prefer N*NB)
559 *
560  irows = ihi + 1 - ilo
561  icols = n + 1 - ilo
562  itau = 1
563  iwrk = itau + irows
564  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
565  $ work( iwrk ), lwork+1-iwrk, ierr )
566 *
567 * Apply the unitary transformation to matrix A
568 * (Complex Workspace: need N, prefer N*NB)
569 *
570  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
571  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
572  $ lwork+1-iwrk, ierr )
573 *
574 * Initialize VSL
575 * (Complex Workspace: need N, prefer N*NB)
576 *
577  IF( ilvsl ) THEN
578  CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
579  IF( irows.GT.1 ) THEN
580  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
581  $ vsl( ilo+1, ilo ), ldvsl )
582  END IF
583  CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
584  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
585  END IF
586 *
587 * Initialize VSR
588 *
589  IF( ilvsr )
590  $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
591 *
592 * Reduce to generalized Hessenberg form
593 * (Workspace: none needed)
594 *
595  CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
596  $ ldvsl, vsr, ldvsr, ierr )
597 *
598  sdim = 0
599 *
600 * Perform QZ algorithm, computing Schur vectors if desired
601 * (Complex Workspace: need N)
602 * (Real Workspace: need N)
603 *
604  iwrk = itau
605  CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
606  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
607  $ lwork+1-iwrk, rwork( irwrk ), ierr )
608  IF( ierr.NE.0 ) THEN
609  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
610  info = ierr
611  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
612  info = ierr - n
613  ELSE
614  info = n + 1
615  END IF
616  GO TO 40
617  END IF
618 *
619 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
620 * condition number(s)
621 *
622  IF( wantst ) THEN
623 *
624 * Undo scaling on eigenvalues before SELCTGing
625 *
626  IF( ilascl )
627  $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
628  IF( ilbscl )
629  $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
630 *
631 * Select eigenvalues
632 *
633  DO 10 i = 1, n
634  bwork( i ) = selctg( alpha( i ), beta( i ) )
635  10 CONTINUE
636 *
637 * Reorder eigenvalues, transform Generalized Schur vectors, and
638 * compute reciprocal condition numbers
639 * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
640 * otherwise, need 1 )
641 *
642  CALL ztgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
643  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, sdim, pl, pr,
644  $ dif, work( iwrk ), lwork-iwrk+1, iwork, liwork,
645  $ ierr )
646 *
647  IF( ijob.GE.1 )
648  $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
649  IF( ierr.EQ.-21 ) THEN
650 *
651 * not enough complex workspace
652 *
653  info = -21
654  ELSE
655  IF( ijob.EQ.1 .OR. ijob.EQ.4 ) THEN
656  rconde( 1 ) = pl
657  rconde( 2 ) = pr
658  END IF
659  IF( ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
660  rcondv( 1 ) = dif( 1 )
661  rcondv( 2 ) = dif( 2 )
662  END IF
663  IF( ierr.EQ.1 )
664  $ info = n + 3
665  END IF
666 *
667  END IF
668 *
669 * Apply permutation to VSL and VSR
670 * (Workspace: none needed)
671 *
672  IF( ilvsl )
673  $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
674  $ rwork( iright ), n, vsl, ldvsl, ierr )
675 *
676  IF( ilvsr )
677  $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
678  $ rwork( iright ), n, vsr, ldvsr, ierr )
679 *
680 * Undo scaling
681 *
682  IF( ilascl ) THEN
683  CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
684  CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
685  END IF
686 *
687  IF( ilbscl ) THEN
688  CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
689  CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
690  END IF
691 *
692  IF( wantst ) THEN
693 *
694 * Check if reordering is correct
695 *
696  lastsl = .true.
697  sdim = 0
698  DO 30 i = 1, n
699  cursl = selctg( alpha( i ), beta( i ) )
700  IF( cursl )
701  $ sdim = sdim + 1
702  IF( cursl .AND. .NOT.lastsl )
703  $ info = n + 2
704  lastsl = cursl
705  30 CONTINUE
706 *
707  END IF
708 *
709  40 CONTINUE
710 *
711  work( 1 ) = maxwrk
712  iwork( 1 ) = liwmin
713 *
714  RETURN
715 *
716 * End of ZGGESX
717 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine ztgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
ZTGSEN
Definition: ztgsen.f:435
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zggev ( character  JOBVL,
character  JOBVR,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, and optionally, the left and/or
 right generalized eigenvectors.

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar
 lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
 singular. It is usually represented as the pair (alpha,beta), as
 there is a reasonable interpretation for beta=0, and even for both
 being zero.

 The right generalized eigenvector v(j) corresponding to the
 generalized eigenvalue lambda(j) of (A,B) satisfies

              A * v(j) = lambda(j) * B * v(j).

 The left generalized eigenvector u(j) corresponding to the
 generalized eigenvalues lambda(j) of (A,B) satisfies

              u(j)**H * A = lambda(j) * u(j)**H * B

 where u(j)**H is the conjugate-transpose of u(j).
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors.
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the matrix A in the pair (A,B).
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the matrix B in the pair (A,B).
          On exit, B has been overwritten.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left generalized eigenvectors u(j) are
          stored one after another in the columns of VL, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right generalized eigenvectors v(j) are
          stored one after another in the columns of VR, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[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).
          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (8*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          =1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHA(j) and BETA(j) should be
                correct for j=INFO+1,...,N.
          > N:  =N+1: other then QZ iteration failed in DHGEQZ,
                =N+2: error return from DTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 219 of file zggev.f.

219 *
220 * -- LAPACK driver routine (version 3.4.1) --
221 * -- LAPACK is a software package provided by Univ. of Tennessee, --
222 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223 * April 2012
224 *
225 * .. Scalar Arguments ..
226  CHARACTER jobvl, jobvr
227  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
228 * ..
229 * .. Array Arguments ..
230  DOUBLE PRECISION rwork( * )
231  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
232  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
233  $ work( * )
234 * ..
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239  DOUBLE PRECISION zero, one
240  parameter( zero = 0.0d0, one = 1.0d0 )
241  COMPLEX*16 czero, cone
242  parameter( czero = ( 0.0d0, 0.0d0 ),
243  $ cone = ( 1.0d0, 0.0d0 ) )
244 * ..
245 * .. Local Scalars ..
246  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery
247  CHARACTER chtemp
248  INTEGER icols, ierr, ihi, ijobvl, ijobvr, ileft, ilo,
249  $ in, iright, irows, irwrk, itau, iwrk, jc, jr,
250  $ lwkmin, lwkopt
251  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
252  $ smlnum, temp
253  COMPLEX*16 x
254 * ..
255 * .. Local Arrays ..
256  LOGICAL ldumma( 1 )
257 * ..
258 * .. External Subroutines ..
259  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghrd,
261  $ zunmqr
262 * ..
263 * .. External Functions ..
264  LOGICAL lsame
265  INTEGER ilaenv
266  DOUBLE PRECISION dlamch, zlange
267  EXTERNAL lsame, ilaenv, dlamch, zlange
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC abs, dble, dimag, max, sqrt
271 * ..
272 * .. Statement Functions ..
273  DOUBLE PRECISION abs1
274 * ..
275 * .. Statement Function definitions ..
276  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
277 * ..
278 * .. Executable Statements ..
279 *
280 * Decode the input arguments
281 *
282  IF( lsame( jobvl, 'N' ) ) THEN
283  ijobvl = 1
284  ilvl = .false.
285  ELSE IF( lsame( jobvl, 'V' ) ) THEN
286  ijobvl = 2
287  ilvl = .true.
288  ELSE
289  ijobvl = -1
290  ilvl = .false.
291  END IF
292 *
293  IF( lsame( jobvr, 'N' ) ) THEN
294  ijobvr = 1
295  ilvr = .false.
296  ELSE IF( lsame( jobvr, 'V' ) ) THEN
297  ijobvr = 2
298  ilvr = .true.
299  ELSE
300  ijobvr = -1
301  ilvr = .false.
302  END IF
303  ilv = ilvl .OR. ilvr
304 *
305 * Test the input arguments
306 *
307  info = 0
308  lquery = ( lwork.EQ.-1 )
309  IF( ijobvl.LE.0 ) THEN
310  info = -1
311  ELSE IF( ijobvr.LE.0 ) THEN
312  info = -2
313  ELSE IF( n.LT.0 ) THEN
314  info = -3
315  ELSE IF( lda.LT.max( 1, n ) ) THEN
316  info = -5
317  ELSE IF( ldb.LT.max( 1, n ) ) THEN
318  info = -7
319  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
320  info = -11
321  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
322  info = -13
323  END IF
324 *
325 * Compute workspace
326 * (Note: Comments in the code beginning "Workspace:" describe the
327 * minimal amount of workspace needed at that point in the code,
328 * as well as the preferred amount for good performance.
329 * NB refers to the optimal block size for the immediately
330 * following subroutine, as returned by ILAENV. The workspace is
331 * computed assuming ILO = 1 and IHI = N, the worst case.)
332 *
333  IF( info.EQ.0 ) THEN
334  lwkmin = max( 1, 2*n )
335  lwkopt = max( 1, n + n*ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
336  lwkopt = max( lwkopt, n +
337  $ n*ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, 0 ) )
338  IF( ilvl ) THEN
339  lwkopt = max( lwkopt, n +
340  $ n*ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, -1 ) )
341  END IF
342  work( 1 ) = lwkopt
343 *
344  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
345  $ info = -15
346  END IF
347 *
348  IF( info.NE.0 ) THEN
349  CALL xerbla( 'ZGGEV ', -info )
350  RETURN
351  ELSE IF( lquery ) THEN
352  RETURN
353  END IF
354 *
355 * Quick return if possible
356 *
357  IF( n.EQ.0 )
358  $ RETURN
359 *
360 * Get machine constants
361 *
362  eps = dlamch( 'E' )*dlamch( 'B' )
363  smlnum = dlamch( 'S' )
364  bignum = one / smlnum
365  CALL dlabad( smlnum, bignum )
366  smlnum = sqrt( smlnum ) / eps
367  bignum = one / smlnum
368 *
369 * Scale A if max element outside range [SMLNUM,BIGNUM]
370 *
371  anrm = zlange( 'M', n, n, a, lda, rwork )
372  ilascl = .false.
373  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
374  anrmto = smlnum
375  ilascl = .true.
376  ELSE IF( anrm.GT.bignum ) THEN
377  anrmto = bignum
378  ilascl = .true.
379  END IF
380  IF( ilascl )
381  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
382 *
383 * Scale B if max element outside range [SMLNUM,BIGNUM]
384 *
385  bnrm = zlange( 'M', n, n, b, ldb, rwork )
386  ilbscl = .false.
387  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
388  bnrmto = smlnum
389  ilbscl = .true.
390  ELSE IF( bnrm.GT.bignum ) THEN
391  bnrmto = bignum
392  ilbscl = .true.
393  END IF
394  IF( ilbscl )
395  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
396 *
397 * Permute the matrices A, B to isolate eigenvalues if possible
398 * (Real Workspace: need 6*N)
399 *
400  ileft = 1
401  iright = n + 1
402  irwrk = iright + n
403  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
404  $ rwork( iright ), rwork( irwrk ), ierr )
405 *
406 * Reduce B to triangular form (QR decomposition of B)
407 * (Complex Workspace: need N, prefer N*NB)
408 *
409  irows = ihi + 1 - ilo
410  IF( ilv ) THEN
411  icols = n + 1 - ilo
412  ELSE
413  icols = irows
414  END IF
415  itau = 1
416  iwrk = itau + irows
417  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
418  $ work( iwrk ), lwork+1-iwrk, ierr )
419 *
420 * Apply the orthogonal transformation to matrix A
421 * (Complex Workspace: need N, prefer N*NB)
422 *
423  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
424  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
425  $ lwork+1-iwrk, ierr )
426 *
427 * Initialize VL
428 * (Complex Workspace: need N, prefer N*NB)
429 *
430  IF( ilvl ) THEN
431  CALL zlaset( 'Full', n, n, czero, cone, vl, ldvl )
432  IF( irows.GT.1 ) THEN
433  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
434  $ vl( ilo+1, ilo ), ldvl )
435  END IF
436  CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
437  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
438  END IF
439 *
440 * Initialize VR
441 *
442  IF( ilvr )
443  $ CALL zlaset( 'Full', n, n, czero, cone, vr, ldvr )
444 *
445 * Reduce to generalized Hessenberg form
446 *
447  IF( ilv ) THEN
448 *
449 * Eigenvectors requested -- work on whole matrix.
450 *
451  CALL zgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
452  $ ldvl, vr, ldvr, ierr )
453  ELSE
454  CALL zgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
455  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
456  END IF
457 *
458 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
459 * Schur form and Schur vectors)
460 * (Complex Workspace: need N)
461 * (Real Workspace: need N)
462 *
463  iwrk = itau
464  IF( ilv ) THEN
465  chtemp = 'S'
466  ELSE
467  chtemp = 'E'
468  END IF
469  CALL zhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
470  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
471  $ lwork+1-iwrk, rwork( irwrk ), ierr )
472  IF( ierr.NE.0 ) THEN
473  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
474  info = ierr
475  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
476  info = ierr - n
477  ELSE
478  info = n + 1
479  END IF
480  GO TO 70
481  END IF
482 *
483 * Compute Eigenvectors
484 * (Real Workspace: need 2*N)
485 * (Complex Workspace: need 2*N)
486 *
487  IF( ilv ) THEN
488  IF( ilvl ) THEN
489  IF( ilvr ) THEN
490  chtemp = 'B'
491  ELSE
492  chtemp = 'L'
493  END IF
494  ELSE
495  chtemp = 'R'
496  END IF
497 *
498  CALL ztgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
499  $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
500  $ ierr )
501  IF( ierr.NE.0 ) THEN
502  info = n + 2
503  GO TO 70
504  END IF
505 *
506 * Undo balancing on VL and VR and normalization
507 * (Workspace: none needed)
508 *
509  IF( ilvl ) THEN
510  CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
511  $ rwork( iright ), n, vl, ldvl, ierr )
512  DO 30 jc = 1, n
513  temp = zero
514  DO 10 jr = 1, n
515  temp = max( temp, abs1( vl( jr, jc ) ) )
516  10 CONTINUE
517  IF( temp.LT.smlnum )
518  $ GO TO 30
519  temp = one / temp
520  DO 20 jr = 1, n
521  vl( jr, jc ) = vl( jr, jc )*temp
522  20 CONTINUE
523  30 CONTINUE
524  END IF
525  IF( ilvr ) THEN
526  CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
527  $ rwork( iright ), n, vr, ldvr, ierr )
528  DO 60 jc = 1, n
529  temp = zero
530  DO 40 jr = 1, n
531  temp = max( temp, abs1( vr( jr, jc ) ) )
532  40 CONTINUE
533  IF( temp.LT.smlnum )
534  $ GO TO 60
535  temp = one / temp
536  DO 50 jr = 1, n
537  vr( jr, jc ) = vr( jr, jc )*temp
538  50 CONTINUE
539  60 CONTINUE
540  END IF
541  END IF
542 *
543 * Undo scaling if necessary
544 *
545  70 CONTINUE
546 *
547  IF( ilascl )
548  $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
549 *
550  IF( ilbscl )
551  $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
552 *
553  work( 1 ) = lwkopt
554  RETURN
555 *
556 * End of ZGGEV
557 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
Definition: ztgevc.f:221
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zggev3 ( character  JOBVL,
character  JOBVR,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)

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

Purpose:
 ZGGEV3 computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, and optionally, the left and/or
 right generalized eigenvectors.

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar
 lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
 singular. It is usually represented as the pair (alpha,beta), as
 there is a reasonable interpretation for beta=0, and even for both
 being zero.

 The right generalized eigenvector v(j) corresponding to the
 generalized eigenvalue lambda(j) of (A,B) satisfies

              A * v(j) = lambda(j) * B * v(j).

 The left generalized eigenvector u(j) corresponding to the
 generalized eigenvalues lambda(j) of (A,B) satisfies

              u(j)**H * A = lambda(j) * u(j)**H * B

 where u(j)**H is the conjugate-transpose of u(j).
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors.
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the matrix A in the pair (A,B).
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the matrix B in the pair (A,B).
          On exit, B has been overwritten.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left generalized eigenvectors u(j) are
          stored one after another in the columns of VL, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right generalized eigenvectors v(j) are
          stored one after another in the columns of VR, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[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.

          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (8*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          =1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHA(j) and BETA(j) should be
                correct for j=INFO+1,...,N.
          > N:  =N+1: other then QZ iteration failed in DHGEQZ,
                =N+2: error return from DTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
January 2015

Definition at line 218 of file zggev3.f.

218 *
219 * -- LAPACK driver routine (version 3.6.0) --
220 * -- LAPACK is a software package provided by Univ. of Tennessee, --
221 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222 * January 2015
223 *
224 * .. Scalar Arguments ..
225  CHARACTER jobvl, jobvr
226  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
227 * ..
228 * .. Array Arguments ..
229  DOUBLE PRECISION rwork( * )
230  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
231  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
232  $ work( * )
233 * ..
234 *
235 * =====================================================================
236 *
237 * .. Parameters ..
238  DOUBLE PRECISION zero, one
239  parameter( zero = 0.0d0, one = 1.0d0 )
240  COMPLEX*16 czero, cone
241  parameter( czero = ( 0.0d0, 0.0d0 ),
242  $ cone = ( 1.0d0, 0.0d0 ) )
243 * ..
244 * .. Local Scalars ..
245  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery
246  CHARACTER chtemp
247  INTEGER icols, ierr, ihi, ijobvl, ijobvr, ileft, ilo,
248  $ in, iright, irows, irwrk, itau, iwrk, jc, jr,
249  $ lwkopt
250  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
251  $ smlnum, temp
252  COMPLEX*16 x
253 * ..
254 * .. Local Arrays ..
255  LOGICAL ldumma( 1 )
256 * ..
257 * .. External Subroutines ..
258  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghd3,
260  $ zunmqr
261 * ..
262 * .. External Functions ..
263  LOGICAL lsame
264  DOUBLE PRECISION dlamch, zlange
265  EXTERNAL lsame, dlamch, zlange
266 * ..
267 * .. Intrinsic Functions ..
268  INTRINSIC abs, dble, dimag, max, sqrt
269 * ..
270 * .. Statement Functions ..
271  DOUBLE PRECISION abs1
272 * ..
273 * .. Statement Function definitions ..
274  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
275 * ..
276 * .. Executable Statements ..
277 *
278 * Decode the input arguments
279 *
280  IF( lsame( jobvl, 'N' ) ) THEN
281  ijobvl = 1
282  ilvl = .false.
283  ELSE IF( lsame( jobvl, 'V' ) ) THEN
284  ijobvl = 2
285  ilvl = .true.
286  ELSE
287  ijobvl = -1
288  ilvl = .false.
289  END IF
290 *
291  IF( lsame( jobvr, 'N' ) ) THEN
292  ijobvr = 1
293  ilvr = .false.
294  ELSE IF( lsame( jobvr, 'V' ) ) THEN
295  ijobvr = 2
296  ilvr = .true.
297  ELSE
298  ijobvr = -1
299  ilvr = .false.
300  END IF
301  ilv = ilvl .OR. ilvr
302 *
303 * Test the input arguments
304 *
305  info = 0
306  lquery = ( lwork.EQ.-1 )
307  IF( ijobvl.LE.0 ) THEN
308  info = -1
309  ELSE IF( ijobvr.LE.0 ) THEN
310  info = -2
311  ELSE IF( n.LT.0 ) THEN
312  info = -3
313  ELSE IF( lda.LT.max( 1, n ) ) THEN
314  info = -5
315  ELSE IF( ldb.LT.max( 1, n ) ) THEN
316  info = -7
317  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
318  info = -11
319  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
320  info = -13
321  ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
322  info = -15
323  END IF
324 *
325 * Compute workspace
326 *
327  IF( info.EQ.0 ) THEN
328  CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
329  lwkopt = max( 1, n+int( work( 1 ) ) )
330  CALL zunmqr( 'L', 'C', n, n, n, b, ldb, work, a, lda, work,
331  $ -1, ierr )
332  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
333  IF( ilvl ) THEN
334  CALL zungqr( n, n, n, vl, ldvl, work, work, -1, ierr )
335  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
336  END IF
337  IF( ilv ) THEN
338  CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
339  $ ldvl, vr, ldvr, work, -1, ierr )
340  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
341  CALL zhgeqz( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
342  $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
343  $ work, ierr )
344  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
345  ELSE
346  CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
347  $ ldvl, vr, ldvr, work, -1, ierr )
348  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
349  CALL zhgeqz( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
350  $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
351  $ work, ierr )
352  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
353  END IF
354  work( 1 ) = dcmplx( lwkopt )
355  END IF
356 *
357  IF( info.NE.0 ) THEN
358  CALL xerbla( 'ZGGEV3 ', -info )
359  RETURN
360  ELSE IF( lquery ) THEN
361  RETURN
362  END IF
363 *
364 * Quick return if possible
365 *
366  IF( n.EQ.0 )
367  $ RETURN
368 *
369 * Get machine constants
370 *
371  eps = dlamch( 'E' )*dlamch( 'B' )
372  smlnum = dlamch( 'S' )
373  bignum = one / smlnum
374  CALL dlabad( smlnum, bignum )
375  smlnum = sqrt( smlnum ) / eps
376  bignum = one / smlnum
377 *
378 * Scale A if max element outside range [SMLNUM,BIGNUM]
379 *
380  anrm = zlange( 'M', n, n, a, lda, rwork )
381  ilascl = .false.
382  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
383  anrmto = smlnum
384  ilascl = .true.
385  ELSE IF( anrm.GT.bignum ) THEN
386  anrmto = bignum
387  ilascl = .true.
388  END IF
389  IF( ilascl )
390  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
391 *
392 * Scale B if max element outside range [SMLNUM,BIGNUM]
393 *
394  bnrm = zlange( 'M', n, n, b, ldb, rwork )
395  ilbscl = .false.
396  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
397  bnrmto = smlnum
398  ilbscl = .true.
399  ELSE IF( bnrm.GT.bignum ) THEN
400  bnrmto = bignum
401  ilbscl = .true.
402  END IF
403  IF( ilbscl )
404  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
405 *
406 * Permute the matrices A, B to isolate eigenvalues if possible
407 *
408  ileft = 1
409  iright = n + 1
410  irwrk = iright + n
411  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
412  $ rwork( iright ), rwork( irwrk ), ierr )
413 *
414 * Reduce B to triangular form (QR decomposition of B)
415 *
416  irows = ihi + 1 - ilo
417  IF( ilv ) THEN
418  icols = n + 1 - ilo
419  ELSE
420  icols = irows
421  END IF
422  itau = 1
423  iwrk = itau + irows
424  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
425  $ work( iwrk ), lwork+1-iwrk, ierr )
426 *
427 * Apply the orthogonal transformation to matrix A
428 *
429  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
430  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
431  $ lwork+1-iwrk, ierr )
432 *
433 * Initialize VL
434 *
435  IF( ilvl ) THEN
436  CALL zlaset( 'Full', n, n, czero, cone, vl, ldvl )
437  IF( irows.GT.1 ) THEN
438  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
439  $ vl( ilo+1, ilo ), ldvl )
440  END IF
441  CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
442  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
443  END IF
444 *
445 * Initialize VR
446 *
447  IF( ilvr )
448  $ CALL zlaset( 'Full', n, n, czero, cone, vr, ldvr )
449 *
450 * Reduce to generalized Hessenberg form
451 *
452  IF( ilv ) THEN
453 *
454 * Eigenvectors requested -- work on whole matrix.
455 *
456  CALL zgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
457  $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
458  ELSE
459  CALL zgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
460  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
461  $ work( iwrk ), lwork+1-iwrk, ierr )
462  END IF
463 *
464 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
465 * Schur form and Schur vectors)
466 *
467  iwrk = itau
468  IF( ilv ) THEN
469  chtemp = 'S'
470  ELSE
471  chtemp = 'E'
472  END IF
473  CALL zhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
474  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
475  $ lwork+1-iwrk, rwork( irwrk ), ierr )
476  IF( ierr.NE.0 ) THEN
477  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
478  info = ierr
479  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
480  info = ierr - n
481  ELSE
482  info = n + 1
483  END IF
484  GO TO 70
485  END IF
486 *
487 * Compute Eigenvectors
488 *
489  IF( ilv ) THEN
490  IF( ilvl ) THEN
491  IF( ilvr ) THEN
492  chtemp = 'B'
493  ELSE
494  chtemp = 'L'
495  END IF
496  ELSE
497  chtemp = 'R'
498  END IF
499 *
500  CALL ztgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
501  $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
502  $ ierr )
503  IF( ierr.NE.0 ) THEN
504  info = n + 2
505  GO TO 70
506  END IF
507 *
508 * Undo balancing on VL and VR and normalization
509 *
510  IF( ilvl ) THEN
511  CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
512  $ rwork( iright ), n, vl, ldvl, ierr )
513  DO 30 jc = 1, n
514  temp = zero
515  DO 10 jr = 1, n
516  temp = max( temp, abs1( vl( jr, jc ) ) )
517  10 CONTINUE
518  IF( temp.LT.smlnum )
519  $ GO TO 30
520  temp = one / temp
521  DO 20 jr = 1, n
522  vl( jr, jc ) = vl( jr, jc )*temp
523  20 CONTINUE
524  30 CONTINUE
525  END IF
526  IF( ilvr ) THEN
527  CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
528  $ rwork( iright ), n, vr, ldvr, ierr )
529  DO 60 jc = 1, n
530  temp = zero
531  DO 40 jr = 1, n
532  temp = max( temp, abs1( vr( jr, jc ) ) )
533  40 CONTINUE
534  IF( temp.LT.smlnum )
535  $ GO TO 60
536  temp = one / temp
537  DO 50 jr = 1, n
538  vr( jr, jc ) = vr( jr, jc )*temp
539  50 CONTINUE
540  60 CONTINUE
541  END IF
542  END IF
543 *
544 * Undo scaling if necessary
545 *
546  70 CONTINUE
547 *
548  IF( ilascl )
549  $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
550 *
551  IF( ilbscl )
552  $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
553 *
554  work( 1 ) = dcmplx( lwkopt )
555  RETURN
556 *
557 * End of ZGGEV3
558 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
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:141
subroutine zgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
ZGGHD3
Definition: zgghd3.f:229
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
Definition: ztgevc.f:221
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zggevx ( character  BALANC,
character  JOBVL,
character  JOBVR,
character  SENSE,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
integer  ILO,
integer  IHI,
double precision, dimension( * )  LSCALE,
double precision, dimension( * )  RSCALE,
double precision  ABNRM,
double precision  BBNRM,
double precision, dimension( * )  RCONDE,
double precision, dimension( * )  RCONDV,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

ZGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B) the generalized eigenvalues, and optionally, the left and/or
 right generalized eigenvectors.

 Optionally, it also computes a balancing transformation to improve
 the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
 LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
 the eigenvalues (RCONDE), and reciprocal condition numbers for the
 right eigenvectors (RCONDV).

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar
 lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
 singular. It is usually represented as the pair (alpha,beta), as
 there is a reasonable interpretation for beta=0, and even for both
 being zero.

 The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
 of (A,B) satisfies
                  A * v(j) = lambda(j) * B * v(j) .
 The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
 of (A,B) satisfies
                  u(j)**H * A  = lambda(j) * u(j)**H * B.
 where u(j)**H is the conjugate-transpose of u(j).
Parameters
[in]BALANC
          BALANC is CHARACTER*1
          Specifies the balance option to be performed:
          = 'N':  do not diagonally scale or permute;
          = 'P':  permute only;
          = 'S':  scale only;
          = 'B':  both permute and scale.
          Computed reciprocal condition numbers will be for the
          matrices after permuting and/or balancing. Permuting does
          not change condition numbers (in exact arithmetic), but
          balancing does.
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors.
[in]SENSE
          SENSE is CHARACTER*1
          Determines which reciprocal condition numbers are computed.
          = 'N': none are computed;
          = 'E': computed for eigenvalues only;
          = 'V': computed for eigenvectors only;
          = 'B': computed for eigenvalues and eigenvectors.
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the matrix A in the pair (A,B).
          On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
          or both, then A contains the first part of the complex Schur
          form of the "balanced" versions of the input A and B.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the matrix B in the pair (A,B).
          On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
          or both, then B contains the second part of the complex
          Schur form of the "balanced" versions of the input A and B.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
          eigenvalues.

          Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio ALPHA/BETA.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left generalized eigenvectors u(j) are
          stored one after another in the columns of VL, in the same
          order as their eigenvalues.
          Each eigenvector will be scaled so the largest component
          will have abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right generalized eigenvectors v(j) are
          stored one after another in the columns of VR, in the same
          order as their eigenvalues.
          Each eigenvector will be scaled so the largest component
          will have abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[out]ILO
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are integer values such that on exit
          A(i,j) = 0 and B(i,j) = 0 if i > j and
          j = 1,...,ILO-1 or i = IHI+1,...,N.
          If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
[out]LSCALE
          LSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the left side of A and B.  If PL(j) is the index of the
          row interchanged with row j, and DL(j) is the scaling
          factor applied to row j, then
            LSCALE(j) = PL(j)  for j = 1,...,ILO-1
                      = DL(j)  for j = ILO,...,IHI
                      = PL(j)  for j = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]RSCALE
          RSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the right side of A and B.  If PR(j) is the index of the
          column interchanged with column j, and DR(j) is the scaling
          factor applied to column j, then
            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
                      = DR(j)  for j = ILO,...,IHI
                      = PR(j)  for j = IHI+1,...,N
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]ABNRM
          ABNRM is DOUBLE PRECISION
          The one-norm of the balanced matrix A.
[out]BBNRM
          BBNRM is DOUBLE PRECISION
          The one-norm of the balanced matrix B.
[out]RCONDE
          RCONDE is DOUBLE PRECISION array, dimension (N)
          If SENSE = 'E' or 'B', the reciprocal condition numbers of
          the eigenvalues, stored in consecutive elements of the array.
          If SENSE = 'N' or 'V', RCONDE is not referenced.
[out]RCONDV
          RCONDV is DOUBLE PRECISION array, dimension (N)
          If JOB = 'V' or 'B', the estimated reciprocal condition
          numbers of the eigenvectors, stored in consecutive elements
          of the array. If the eigenvalues cannot be reordered to
          compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
          when the true value would be very small anyway.
          If SENSE = 'N' or 'E', RCONDV is not referenced.
[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).
          If SENSE = 'E', LWORK >= max(1,4*N).
          If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).

          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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (lrwork)
          lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
          and at least max(1,2*N) otherwise.
          Real workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (N+2)
          If SENSE = 'E', IWORK is not referenced.
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
          If SENSE = 'N', BWORK is not referenced.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          = 1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHA(j) and BETA(j) should be correct
                for j=INFO+1,...,N.
          > N:  =N+1: other than QZ iteration failed in ZHGEQZ.
                =N+2: error return from ZTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012
Further Details:
  Balancing a matrix pair (A,B) includes, first, permuting rows and
  columns to isolate eigenvalues, second, applying diagonal similarity
  transformation to the rows and columns to make the rows and columns
  as close in norm as possible. The computed reciprocal condition
  numbers correspond to the balanced matrix. Permuting rows and columns
  will not change the condition numbers (in exact arithmetic) but
  diagonal scaling will.  For further explanation of balancing, see
  section 4.11.1.2 of LAPACK Users' Guide.

  An approximate error bound on the chordal distance between the i-th
  computed generalized eigenvalue w and the corresponding exact
  eigenvalue lambda is

       chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)

  An approximate error bound for the angle between the i-th computed
  eigenvector VL(i) or VR(i) is given by

       EPS * norm(ABNRM, BBNRM) / DIF(i).

  For further explanation of the reciprocal condition numbers RCONDE
  and RCONDV, see section 4.11 of LAPACK User's Guide.

Definition at line 376 of file zggevx.f.

376 *
377 * -- LAPACK driver routine (version 3.4.1) --
378 * -- LAPACK is a software package provided by Univ. of Tennessee, --
379 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
380 * April 2012
381 *
382 * .. Scalar Arguments ..
383  CHARACTER balanc, jobvl, jobvr, sense
384  INTEGER ihi, ilo, info, lda, ldb, ldvl, ldvr, lwork, n
385  DOUBLE PRECISION abnrm, bbnrm
386 * ..
387 * .. Array Arguments ..
388  LOGICAL bwork( * )
389  INTEGER iwork( * )
390  DOUBLE PRECISION lscale( * ), rconde( * ), rcondv( * ),
391  $ rscale( * ), rwork( * )
392  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
393  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
394  $ work( * )
395 * ..
396 *
397 * =====================================================================
398 *
399 * .. Parameters ..
400  DOUBLE PRECISION zero, one
401  parameter( zero = 0.0d+0, one = 1.0d+0 )
402  COMPLEX*16 czero, cone
403  parameter( czero = ( 0.0d+0, 0.0d+0 ),
404  $ cone = ( 1.0d+0, 0.0d+0 ) )
405 * ..
406 * .. Local Scalars ..
407  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery, noscl,
408  $ wantsb, wantse, wantsn, wantsv
409  CHARACTER chtemp
410  INTEGER i, icols, ierr, ijobvl, ijobvr, in, irows,
411  $ itau, iwrk, iwrk1, j, jc, jr, m, maxwrk, minwrk
412  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
413  $ smlnum, temp
414  COMPLEX*16 x
415 * ..
416 * .. Local Arrays ..
417  LOGICAL ldumma( 1 )
418 * ..
419 * .. External Subroutines ..
420  EXTERNAL dlabad, dlascl, xerbla, zgeqrf, zggbak, zggbal,
422  $ ztgsna, zungqr, zunmqr
423 * ..
424 * .. External Functions ..
425  LOGICAL lsame
426  INTEGER ilaenv
427  DOUBLE PRECISION dlamch, zlange
428  EXTERNAL lsame, ilaenv, dlamch, zlange
429 * ..
430 * .. Intrinsic Functions ..
431  INTRINSIC abs, dble, dimag, max, sqrt
432 * ..
433 * .. Statement Functions ..
434  DOUBLE PRECISION abs1
435 * ..
436 * .. Statement Function definitions ..
437  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
438 * ..
439 * .. Executable Statements ..
440 *
441 * Decode the input arguments
442 *
443  IF( lsame( jobvl, 'N' ) ) THEN
444  ijobvl = 1
445  ilvl = .false.
446  ELSE IF( lsame( jobvl, 'V' ) ) THEN
447  ijobvl = 2
448  ilvl = .true.
449  ELSE
450  ijobvl = -1
451  ilvl = .false.
452  END IF
453 *
454  IF( lsame( jobvr, 'N' ) ) THEN
455  ijobvr = 1
456  ilvr = .false.
457  ELSE IF( lsame( jobvr, 'V' ) ) THEN
458  ijobvr = 2
459  ilvr = .true.
460  ELSE
461  ijobvr = -1
462  ilvr = .false.
463  END IF
464  ilv = ilvl .OR. ilvr
465 *
466  noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
467  wantsn = lsame( sense, 'N' )
468  wantse = lsame( sense, 'E' )
469  wantsv = lsame( sense, 'V' )
470  wantsb = lsame( sense, 'B' )
471 *
472 * Test the input arguments
473 *
474  info = 0
475  lquery = ( lwork.EQ.-1 )
476  IF( .NOT.( noscl .OR. lsame( balanc,'S' ) .OR.
477  $ lsame( balanc, 'B' ) ) ) THEN
478  info = -1
479  ELSE IF( ijobvl.LE.0 ) THEN
480  info = -2
481  ELSE IF( ijobvr.LE.0 ) THEN
482  info = -3
483  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
484  $ THEN
485  info = -4
486  ELSE IF( n.LT.0 ) THEN
487  info = -5
488  ELSE IF( lda.LT.max( 1, n ) ) THEN
489  info = -7
490  ELSE IF( ldb.LT.max( 1, n ) ) THEN
491  info = -9
492  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
493  info = -13
494  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
495  info = -15
496  END IF
497 *
498 * Compute workspace
499 * (Note: Comments in the code beginning "Workspace:" describe the
500 * minimal amount of workspace needed at that point in the code,
501 * as well as the preferred amount for good performance.
502 * NB refers to the optimal block size for the immediately
503 * following subroutine, as returned by ILAENV. The workspace is
504 * computed assuming ILO = 1 and IHI = N, the worst case.)
505 *
506  IF( info.EQ.0 ) THEN
507  IF( n.EQ.0 ) THEN
508  minwrk = 1
509  maxwrk = 1
510  ELSE
511  minwrk = 2*n
512  IF( wantse ) THEN
513  minwrk = 4*n
514  ELSE IF( wantsv .OR. wantsb ) THEN
515  minwrk = 2*n*( n + 1)
516  END IF
517  maxwrk = minwrk
518  maxwrk = max( maxwrk,
519  $ n + n*ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
520  maxwrk = max( maxwrk,
521  $ n + n*ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, 0 ) )
522  IF( ilvl ) THEN
523  maxwrk = max( maxwrk, n +
524  $ n*ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, 0 ) )
525  END IF
526  END IF
527  work( 1 ) = maxwrk
528 *
529  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
530  info = -25
531  END IF
532  END IF
533 *
534  IF( info.NE.0 ) THEN
535  CALL xerbla( 'ZGGEVX', -info )
536  RETURN
537  ELSE IF( lquery ) THEN
538  RETURN
539  END IF
540 *
541 * Quick return if possible
542 *
543  IF( n.EQ.0 )
544  $ RETURN
545 *
546 * Get machine constants
547 *
548  eps = dlamch( 'P' )
549  smlnum = dlamch( 'S' )
550  bignum = one / smlnum
551  CALL dlabad( smlnum, bignum )
552  smlnum = sqrt( smlnum ) / eps
553  bignum = one / smlnum
554 *
555 * Scale A if max element outside range [SMLNUM,BIGNUM]
556 *
557  anrm = zlange( 'M', n, n, a, lda, rwork )
558  ilascl = .false.
559  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
560  anrmto = smlnum
561  ilascl = .true.
562  ELSE IF( anrm.GT.bignum ) THEN
563  anrmto = bignum
564  ilascl = .true.
565  END IF
566  IF( ilascl )
567  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
568 *
569 * Scale B if max element outside range [SMLNUM,BIGNUM]
570 *
571  bnrm = zlange( 'M', n, n, b, ldb, rwork )
572  ilbscl = .false.
573  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
574  bnrmto = smlnum
575  ilbscl = .true.
576  ELSE IF( bnrm.GT.bignum ) THEN
577  bnrmto = bignum
578  ilbscl = .true.
579  END IF
580  IF( ilbscl )
581  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
582 *
583 * Permute and/or balance the matrix pair (A,B)
584 * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
585 *
586  CALL zggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
587  $ rwork, ierr )
588 *
589 * Compute ABNRM and BBNRM
590 *
591  abnrm = zlange( '1', n, n, a, lda, rwork( 1 ) )
592  IF( ilascl ) THEN
593  rwork( 1 ) = abnrm
594  CALL dlascl( 'G', 0, 0, anrmto, anrm, 1, 1, rwork( 1 ), 1,
595  $ ierr )
596  abnrm = rwork( 1 )
597  END IF
598 *
599  bbnrm = zlange( '1', n, n, b, ldb, rwork( 1 ) )
600  IF( ilbscl ) THEN
601  rwork( 1 ) = bbnrm
602  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, rwork( 1 ), 1,
603  $ ierr )
604  bbnrm = rwork( 1 )
605  END IF
606 *
607 * Reduce B to triangular form (QR decomposition of B)
608 * (Complex Workspace: need N, prefer N*NB )
609 *
610  irows = ihi + 1 - ilo
611  IF( ilv .OR. .NOT.wantsn ) THEN
612  icols = n + 1 - ilo
613  ELSE
614  icols = irows
615  END IF
616  itau = 1
617  iwrk = itau + irows
618  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
619  $ work( iwrk ), lwork+1-iwrk, ierr )
620 *
621 * Apply the unitary transformation to A
622 * (Complex Workspace: need N, prefer N*NB)
623 *
624  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
625  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
626  $ lwork+1-iwrk, ierr )
627 *
628 * Initialize VL and/or VR
629 * (Workspace: need N, prefer N*NB)
630 *
631  IF( ilvl ) THEN
632  CALL zlaset( 'Full', n, n, czero, cone, vl, ldvl )
633  IF( irows.GT.1 ) THEN
634  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
635  $ vl( ilo+1, ilo ), ldvl )
636  END IF
637  CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
638  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
639  END IF
640 *
641  IF( ilvr )
642  $ CALL zlaset( 'Full', n, n, czero, cone, vr, ldvr )
643 *
644 * Reduce to generalized Hessenberg form
645 * (Workspace: none needed)
646 *
647  IF( ilv .OR. .NOT.wantsn ) THEN
648 *
649 * Eigenvectors requested -- work on whole matrix.
650 *
651  CALL zgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
652  $ ldvl, vr, ldvr, ierr )
653  ELSE
654  CALL zgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
655  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
656  END IF
657 *
658 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
659 * Schur forms and Schur vectors)
660 * (Complex Workspace: need N)
661 * (Real Workspace: need N)
662 *
663  iwrk = itau
664  IF( ilv .OR. .NOT.wantsn ) THEN
665  chtemp = 'S'
666  ELSE
667  chtemp = 'E'
668  END IF
669 *
670  CALL zhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
671  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
672  $ lwork+1-iwrk, rwork, ierr )
673  IF( ierr.NE.0 ) THEN
674  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
675  info = ierr
676  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
677  info = ierr - n
678  ELSE
679  info = n + 1
680  END IF
681  GO TO 90
682  END IF
683 *
684 * Compute Eigenvectors and estimate condition numbers if desired
685 * ZTGEVC: (Complex Workspace: need 2*N )
686 * (Real Workspace: need 2*N )
687 * ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
688 * (Integer Workspace: need N+2 )
689 *
690  IF( ilv .OR. .NOT.wantsn ) THEN
691  IF( ilv ) THEN
692  IF( ilvl ) THEN
693  IF( ilvr ) THEN
694  chtemp = 'B'
695  ELSE
696  chtemp = 'L'
697  END IF
698  ELSE
699  chtemp = 'R'
700  END IF
701 *
702  CALL ztgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
703  $ ldvl, vr, ldvr, n, in, work( iwrk ), rwork,
704  $ ierr )
705  IF( ierr.NE.0 ) THEN
706  info = n + 2
707  GO TO 90
708  END IF
709  END IF
710 *
711  IF( .NOT.wantsn ) THEN
712 *
713 * compute eigenvectors (DTGEVC) and estimate condition
714 * numbers (DTGSNA). Note that the definition of the condition
715 * number is not invariant under transformation (u,v) to
716 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
717 * Schur form (S,T), Q and Z are orthogonal matrices. In order
718 * to avoid using extra 2*N*N workspace, we have to
719 * re-calculate eigenvectors and estimate the condition numbers
720 * one at a time.
721 *
722  DO 20 i = 1, n
723 *
724  DO 10 j = 1, n
725  bwork( j ) = .false.
726  10 CONTINUE
727  bwork( i ) = .true.
728 *
729  iwrk = n + 1
730  iwrk1 = iwrk + n
731 *
732  IF( wantse .OR. wantsb ) THEN
733  CALL ztgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
734  $ work( 1 ), n, work( iwrk ), n, 1, m,
735  $ work( iwrk1 ), rwork, ierr )
736  IF( ierr.NE.0 ) THEN
737  info = n + 2
738  GO TO 90
739  END IF
740  END IF
741 *
742  CALL ztgsna( sense, 'S', bwork, n, a, lda, b, ldb,
743  $ work( 1 ), n, work( iwrk ), n, rconde( i ),
744  $ rcondv( i ), 1, m, work( iwrk1 ),
745  $ lwork-iwrk1+1, iwork, ierr )
746 *
747  20 CONTINUE
748  END IF
749  END IF
750 *
751 * Undo balancing on VL and VR and normalization
752 * (Workspace: none needed)
753 *
754  IF( ilvl ) THEN
755  CALL zggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
756  $ ldvl, ierr )
757 *
758  DO 50 jc = 1, n
759  temp = zero
760  DO 30 jr = 1, n
761  temp = max( temp, abs1( vl( jr, jc ) ) )
762  30 CONTINUE
763  IF( temp.LT.smlnum )
764  $ GO TO 50
765  temp = one / temp
766  DO 40 jr = 1, n
767  vl( jr, jc ) = vl( jr, jc )*temp
768  40 CONTINUE
769  50 CONTINUE
770  END IF
771 *
772  IF( ilvr ) THEN
773  CALL zggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
774  $ ldvr, ierr )
775  DO 80 jc = 1, n
776  temp = zero
777  DO 60 jr = 1, n
778  temp = max( temp, abs1( vr( jr, jc ) ) )
779  60 CONTINUE
780  IF( temp.LT.smlnum )
781  $ GO TO 80
782  temp = one / temp
783  DO 70 jr = 1, n
784  vr( jr, jc ) = vr( jr, jc )*temp
785  70 CONTINUE
786  80 CONTINUE
787  END IF
788 *
789 * Undo scaling if necessary
790 *
791  90 CONTINUE
792 *
793  IF( ilascl )
794  $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
795 *
796  IF( ilbscl )
797  $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
798 *
799  work( 1 ) = maxwrk
800  RETURN
801 *
802 * End of ZGGEVX
803 *
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine ztgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
ZTGSNA
Definition: ztgsna.f:313
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
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:141
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
Definition: ztgevc.f:221
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:117
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:141
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286

Here is the call graph for this function:

Here is the caller graph for this function: