OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
complex16
LAPACK » General Matrices » Eigenvalue

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
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
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
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
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
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
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
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)
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
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
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)
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

Detailed Description

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

Function Documentation

◆ zgees()

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.

Definition at line 195 of file zgees.f.

197*
198* -- LAPACK driver routine --
199* -- LAPACK is a software package provided by Univ. of Tennessee, --
200* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201*
202* .. Scalar Arguments ..
203 CHARACTER JOBVS, SORT
204 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
205* ..
206* .. Array Arguments ..
207 LOGICAL BWORK( * )
208 DOUBLE PRECISION RWORK( * )
209 COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
210* ..
211* .. Function Arguments ..
212 LOGICAL SELECT
213 EXTERNAL SELECT
214* ..
215*
216* =====================================================================
217*
218* .. Parameters ..
219 DOUBLE PRECISION ZERO, ONE
220 parameter( zero = 0.0d0, one = 1.0d0 )
221* ..
222* .. Local Scalars ..
223 LOGICAL LQUERY, SCALEA, WANTST, WANTVS
224 INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
225 $ ITAU, IWRK, MAXWRK, MINWRK
226 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
227* ..
228* .. Local Arrays ..
229 DOUBLE PRECISION DUM( 1 )
230* ..
231* .. External Subroutines ..
232 EXTERNAL dlabad, xerbla, zcopy, zgebak, zgebal, zgehrd,
233 $ zhseqr, zlacpy, zlascl, ztrsen, zunghr
234* ..
235* .. External Functions ..
236 LOGICAL LSAME
237 INTEGER ILAENV
238 DOUBLE PRECISION DLAMCH, ZLANGE
239 EXTERNAL lsame, ilaenv, dlamch, zlange
240* ..
241* .. Intrinsic Functions ..
242 INTRINSIC max, sqrt
243* ..
244* .. Executable Statements ..
245*
246* Test the input arguments
247*
248 info = 0
249 lquery = ( lwork.EQ.-1 )
250 wantvs = lsame( jobvs, 'V' )
251 wantst = lsame( sort, 'S' )
252 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
253 info = -1
254 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
255 info = -2
256 ELSE IF( n.LT.0 ) THEN
257 info = -4
258 ELSE IF( lda.LT.max( 1, n ) ) THEN
259 info = -6
260 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
261 info = -10
262 END IF
263*
264* Compute workspace
265* (Note: Comments in the code beginning "Workspace:" describe the
266* minimal amount of workspace needed at that point in the code,
267* as well as the preferred amount for good performance.
268* CWorkspace refers to complex workspace, and RWorkspace to real
269* workspace. NB refers to the optimal block size for the
270* immediately following subroutine, as returned by ILAENV.
271* HSWORK refers to the workspace preferred by ZHSEQR, as
272* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
273* the worst case.)
274*
275 IF( info.EQ.0 ) THEN
276 IF( n.EQ.0 ) THEN
277 minwrk = 1
278 maxwrk = 1
279 ELSE
280 maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
281 minwrk = 2*n
282*
283 CALL zhseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
284 $ work, -1, ieval )
285 hswork = dble( work( 1 ) )
286*
287 IF( .NOT.wantvs ) THEN
288 maxwrk = max( maxwrk, hswork )
289 ELSE
290 maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
291 $ ' ', n, 1, n, -1 ) )
292 maxwrk = max( maxwrk, hswork )
293 END IF
294 END IF
295 work( 1 ) = maxwrk
296*
297 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
298 info = -12
299 END IF
300 END IF
301*
302 IF( info.NE.0 ) THEN
303 CALL xerbla( 'ZGEES ', -info )
304 RETURN
305 ELSE IF( lquery ) THEN
306 RETURN
307 END IF
308*
309* Quick return if possible
310*
311 IF( n.EQ.0 ) THEN
312 sdim = 0
313 RETURN
314 END IF
315*
316* Get machine constants
317*
318 eps = dlamch( 'P' )
319 smlnum = dlamch( 'S' )
320 bignum = one / smlnum
321 CALL dlabad( smlnum, bignum )
322 smlnum = sqrt( smlnum ) / eps
323 bignum = one / smlnum
324*
325* Scale A if max element outside range [SMLNUM,BIGNUM]
326*
327 anrm = zlange( 'M', n, n, a, lda, dum )
328 scalea = .false.
329 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
330 scalea = .true.
331 cscale = smlnum
332 ELSE IF( anrm.GT.bignum ) THEN
333 scalea = .true.
334 cscale = bignum
335 END IF
336 IF( scalea )
337 $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
338*
339* Permute the matrix to make it more nearly triangular
340* (CWorkspace: none)
341* (RWorkspace: need N)
342*
343 ibal = 1
344 CALL zgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
345*
346* Reduce to upper Hessenberg form
347* (CWorkspace: need 2*N, prefer N+N*NB)
348* (RWorkspace: none)
349*
350 itau = 1
351 iwrk = n + itau
352 CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
353 $ lwork-iwrk+1, ierr )
354*
355 IF( wantvs ) THEN
356*
357* Copy Householder vectors to VS
358*
359 CALL zlacpy( 'L', n, n, a, lda, vs, ldvs )
360*
361* Generate unitary matrix in VS
362* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
363* (RWorkspace: none)
364*
365 CALL zunghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
366 $ lwork-iwrk+1, ierr )
367 END IF
368*
369 sdim = 0
370*
371* Perform QR iteration, accumulating Schur vectors in VS if desired
372* (CWorkspace: need 1, prefer HSWORK (see comments) )
373* (RWorkspace: none)
374*
375 iwrk = itau
376 CALL zhseqr( 'S', jobvs, n, ilo, ihi, a, lda, w, vs, ldvs,
377 $ work( iwrk ), lwork-iwrk+1, ieval )
378 IF( ieval.GT.0 )
379 $ info = ieval
380*
381* Sort eigenvalues if desired
382*
383 IF( wantst .AND. info.EQ.0 ) THEN
384 IF( scalea )
385 $ CALL zlascl( 'G', 0, 0, cscale, anrm, n, 1, w, n, ierr )
386 DO 10 i = 1, n
387 bwork( i ) = SELECT( w( i ) )
388 10 CONTINUE
389*
390* Reorder eigenvalues and transform Schur vectors
391* (CWorkspace: none)
392* (RWorkspace: none)
393*
394 CALL ztrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, w, sdim,
395 $ s, sep, work( iwrk ), lwork-iwrk+1, icond )
396 END IF
397*
398 IF( wantvs ) THEN
399*
400* Undo balancing
401* (CWorkspace: none)
402* (RWorkspace: need N)
403*
404 CALL zgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs, ldvs,
405 $ ierr )
406 END IF
407*
408 IF( scalea ) THEN
409*
410* Undo scaling for the Schur form of A
411*
412 CALL zlascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
413 CALL zcopy( n, a, lda+1, w, 1 )
414 END IF
415*
416 work( 1 ) = maxwrk
417 RETURN
418*
419* End of ZGEES
420*
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:115
zgebal
subroutine zgebal(job, n, a, lda, ilo, ihi, scale, info)
ZGEBAL
Definition zgebal.f:162
zgehrd
subroutine zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZGEHRD
Definition zgehrd.f:167
zgebak
subroutine zgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
ZGEBAK
Definition zgebak.f:131
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:143
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zhseqr
subroutine zhseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
ZHSEQR
Definition zhseqr.f:299
ztrsen
subroutine ztrsen(job, compq, select, n, t, ldt, q, ldq, w, m, s, sep, work, lwork, info)
ZTRSEN
Definition ztrsen.f:264
zunghr
subroutine zunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZUNGHR
Definition zunghr.f:126
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21

◆ zgeesx()

subroutine zgeesx ( character jobvs,
character sort,
external select,
character sense,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
integer sdim,
complex*16, dimension( * ) w,
complex*16, dimension( ldvs, * ) vs,
integer ldvs,
double precision rconde,
double precision rcondv,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> ZGEESX computes for an N-by-N complex nonsymmetric matrix A, the
!> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
!> vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
!>
!> Optionally, it also orders the eigenvalues on the diagonal of the
!> Schur form so that selected eigenvalues are at the top left;
!> computes a reciprocal condition number for the average of the
!> selected eigenvalues (RCONDE); and computes a reciprocal condition
!> number for the right invariant subspace corresponding to the
!> selected eigenvalues (RCONDV).  The leading columns of Z form an
!> orthonormal basis for this invariant subspace.
!>
!> For further explanation of the reciprocal condition numbers RCONDE
!> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
!> these quantities are called s and sep respectively).
!>
!> A complex matrix is in Schur form if it is upper triangular.
!>
Parameters
[in]JOBVS
!>          JOBVS is CHARACTER*1
!>          = 'N': Schur vectors are not computed;
!>          = 'V': Schur vectors are computed.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the Schur form.
!>          = 'N': Eigenvalues are not ordered;
!>          = 'S': Eigenvalues are ordered (see SELECT).
!>
[in]SELECT
!>          SELECT is a LOGICAL FUNCTION of one COMPLEX*16 argument
!>          SELECT must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'S', SELECT is used to select eigenvalues to order
!>          to the top left of the Schur form.
!>          If SORT = 'N', SELECT is not referenced.
!>          An eigenvalue W(j) is selected if SELECT(W(j)) is true.
!>
[in]SENSE
!>          SENSE is CHARACTER*1
!>          Determines which reciprocal condition numbers are computed.
!>          = 'N': None are computed;
!>          = 'E': Computed for average of selected eigenvalues only;
!>          = 'V': Computed for selected right invariant subspace only;
!>          = 'B': Computed for both.
!>          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A. N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA, N)
!>          On entry, the N-by-N matrix A.
!>          On exit, A is overwritten by its Schur form T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues for which
!>                         SELECT is true.
!>
[out]W
!>          W is COMPLEX*16 array, dimension (N)
!>          W contains the computed eigenvalues, in the same order
!>          that they appear on the diagonal of the output Schur form T.
!>
[out]VS
!>          VS is COMPLEX*16 array, dimension (LDVS,N)
!>          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
!>          vectors.
!>          If JOBVS = 'N', VS is not referenced.
!>
[in]LDVS
!>          LDVS is INTEGER
!>          The leading dimension of the array VS.  LDVS >= 1, and if
!>          JOBVS = 'V', LDVS >= N.
!>
[out]RCONDE
!>          RCONDE is DOUBLE PRECISION
!>          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
!>          condition number for the average of the selected eigenvalues.
!>          Not referenced if SENSE = 'N' or 'V'.
!>
[out]RCONDV
!>          RCONDV is DOUBLE PRECISION
!>          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
!>          condition number for the selected right invariant subspace.
!>          Not referenced if SENSE = 'N' or 'E'.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,2*N).
!>          Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
!>          where SDIM is the number of selected eigenvalues computed by
!>          this routine.  Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also
!>          that an error is only returned if LWORK < max(1,2*N), but if
!>          SENSE = 'E' or 'V' or 'B' this may not be large enough.
!>          For good performance, LWORK must generally be larger.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates upper bound on the optimal size of the
!>          array WORK, returns this value as the first entry of the WORK
!>          array, and no error message related to LWORK is issued by
!>          XERBLA.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[out]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value.
!>          > 0: if INFO = i, and i is
!>             <= N: the QR algorithm failed to compute all the
!>                   eigenvalues; elements 1:ILO-1 and i+1:N of W
!>                   contain those eigenvalues which have converged; if
!>                   JOBVS = 'V', VS contains the transformation which
!>                   reduces A to its partially converged Schur form.
!>             = N+1: the eigenvalues could not be reordered because some
!>                   eigenvalues were too close to separate (the problem
!>                   is very ill-conditioned);
!>             = N+2: after reordering, roundoff changed values of some
!>                   complex eigenvalues so that leading eigenvalues in
!>                   the Schur form no longer satisfy SELECT=.TRUE.  This
!>                   could also be caused by underflow due to scaling.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 236 of file zgeesx.f.

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

◆ zgeev()

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 i+1:N of W contain eigenvalues which have
!>                converged.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 178 of file zgeev.f.

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

◆ zgeevx()

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.

Definition at line 285 of file zgeevx.f.

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

◆ zgegs()

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.

Definition at line 222 of file zgegs.f.

225*
226* -- LAPACK driver routine --
227* -- LAPACK is a software package provided by Univ. of Tennessee, --
228* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
229*
230* .. Scalar Arguments ..
231 CHARACTER JOBVSL, JOBVSR
232 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
233* ..
234* .. Array Arguments ..
235 DOUBLE PRECISION RWORK( * )
236 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
237 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
238 $ WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 DOUBLE PRECISION ZERO, ONE
245 parameter( zero = 0.0d0, one = 1.0d0 )
246 COMPLEX*16 CZERO, CONE
247 parameter( czero = ( 0.0d0, 0.0d0 ),
248 $ cone = ( 1.0d0, 0.0d0 ) )
249* ..
250* .. Local Scalars ..
251 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
252 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
253 $ IRIGHT, IROWS, IRWORK, ITAU, IWORK, LOPT,
254 $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
255 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
256 $ SAFMIN, SMLNUM
257* ..
258* .. External Subroutines ..
259 EXTERNAL xerbla, zgeqrf, zggbak, zggbal, zgghrd, zhgeqz,
260 $ zlacpy, zlascl, zlaset, zungqr, zunmqr
261* ..
262* .. External Functions ..
263 LOGICAL LSAME
264 INTEGER ILAENV
265 DOUBLE PRECISION DLAMCH, ZLANGE
266 EXTERNAL lsame, ilaenv, dlamch, zlange
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC int, max
270* ..
271* .. Executable Statements ..
272*
273* Decode the input arguments
274*
275 IF( lsame( jobvsl, 'N' ) ) THEN
276 ijobvl = 1
277 ilvsl = .false.
278 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
279 ijobvl = 2
280 ilvsl = .true.
281 ELSE
282 ijobvl = -1
283 ilvsl = .false.
284 END IF
285*
286 IF( lsame( jobvsr, 'N' ) ) THEN
287 ijobvr = 1
288 ilvsr = .false.
289 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
290 ijobvr = 2
291 ilvsr = .true.
292 ELSE
293 ijobvr = -1
294 ilvsr = .false.
295 END IF
296*
297* Test the input arguments
298*
299 lwkmin = max( 2*n, 1 )
300 lwkopt = lwkmin
301 work( 1 ) = lwkopt
302 lquery = ( lwork.EQ.-1 )
303 info = 0
304 IF( ijobvl.LE.0 ) THEN
305 info = -1
306 ELSE IF( ijobvr.LE.0 ) THEN
307 info = -2
308 ELSE IF( n.LT.0 ) THEN
309 info = -3
310 ELSE IF( lda.LT.max( 1, n ) ) THEN
311 info = -5
312 ELSE IF( ldb.LT.max( 1, n ) ) THEN
313 info = -7
314 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
315 info = -11
316 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
317 info = -13
318 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
319 info = -15
320 END IF
321*
322 IF( info.EQ.0 ) THEN
323 nb1 = ilaenv( 1, 'ZGEQRF', ' ', n, n, -1, -1 )
324 nb2 = ilaenv( 1, 'ZUNMQR', ' ', n, n, n, -1 )
325 nb3 = ilaenv( 1, 'ZUNGQR', ' ', n, n, n, -1 )
326 nb = max( nb1, nb2, nb3 )
327 lopt = n*( nb+1 )
328 work( 1 ) = lopt
329 END IF
330*
331 IF( info.NE.0 ) THEN
332 CALL xerbla( 'ZGEGS ', -info )
333 RETURN
334 ELSE IF( lquery ) THEN
335 RETURN
336 END IF
337*
338* Quick return if possible
339*
340 IF( n.EQ.0 )
341 $ RETURN
342*
343* Get machine constants
344*
345 eps = dlamch( 'E' )*dlamch( 'B' )
346 safmin = dlamch( 'S' )
347 smlnum = n*safmin / eps
348 bignum = one / smlnum
349*
350* Scale A if max element outside range [SMLNUM,BIGNUM]
351*
352 anrm = zlange( 'M', n, n, a, lda, rwork )
353 ilascl = .false.
354 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
355 anrmto = smlnum
356 ilascl = .true.
357 ELSE IF( anrm.GT.bignum ) THEN
358 anrmto = bignum
359 ilascl = .true.
360 END IF
361*
362 IF( ilascl ) THEN
363 CALL zlascl( 'G', -1, -1, anrm, anrmto, n, n, a, lda, iinfo )
364 IF( iinfo.NE.0 ) THEN
365 info = n + 9
366 RETURN
367 END IF
368 END IF
369*
370* Scale B if max element outside range [SMLNUM,BIGNUM]
371*
372 bnrm = zlange( 'M', n, n, b, ldb, rwork )
373 ilbscl = .false.
374 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
375 bnrmto = smlnum
376 ilbscl = .true.
377 ELSE IF( bnrm.GT.bignum ) THEN
378 bnrmto = bignum
379 ilbscl = .true.
380 END IF
381*
382 IF( ilbscl ) THEN
383 CALL zlascl( 'G', -1, -1, bnrm, bnrmto, n, n, b, ldb, iinfo )
384 IF( iinfo.NE.0 ) THEN
385 info = n + 9
386 RETURN
387 END IF
388 END IF
389*
390* Permute the matrix to make it more nearly triangular
391*
392 ileft = 1
393 iright = n + 1
394 irwork = iright + n
395 iwork = 1
396 CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
397 $ rwork( iright ), rwork( irwork ), iinfo )
398 IF( iinfo.NE.0 ) THEN
399 info = n + 1
400 GO TO 10
401 END IF
402*
403* Reduce B to triangular form, and initialize VSL and/or VSR
404*
405 irows = ihi + 1 - ilo
406 icols = n + 1 - ilo
407 itau = iwork
408 iwork = itau + irows
409 CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
410 $ work( iwork ), lwork+1-iwork, iinfo )
411 IF( iinfo.GE.0 )
412 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
413 IF( iinfo.NE.0 ) THEN
414 info = n + 2
415 GO TO 10
416 END IF
417*
418 CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
419 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
420 $ lwork+1-iwork, iinfo )
421 IF( iinfo.GE.0 )
422 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
423 IF( iinfo.NE.0 ) THEN
424 info = n + 3
425 GO TO 10
426 END IF
427*
428 IF( ilvsl ) THEN
429 CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
430 CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
431 $ vsl( ilo+1, ilo ), ldvsl )
432 CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
433 $ work( itau ), work( iwork ), lwork+1-iwork,
434 $ iinfo )
435 IF( iinfo.GE.0 )
436 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
437 IF( iinfo.NE.0 ) THEN
438 info = n + 4
439 GO TO 10
440 END IF
441 END IF
442*
443 IF( ilvsr )
444 $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
445*
446* Reduce to generalized Hessenberg form
447*
448 CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
449 $ ldvsl, vsr, ldvsr, iinfo )
450 IF( iinfo.NE.0 ) THEN
451 info = n + 5
452 GO TO 10
453 END IF
454*
455* Perform QZ algorithm, computing Schur vectors if desired
456*
457 iwork = itau
458 CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
459 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwork ),
460 $ lwork+1-iwork, rwork( irwork ), iinfo )
461 IF( iinfo.GE.0 )
462 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
463 IF( iinfo.NE.0 ) THEN
464 IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
465 info = iinfo
466 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
467 info = iinfo - n
468 ELSE
469 info = n + 6
470 END IF
471 GO TO 10
472 END IF
473*
474* Apply permutation to VSL and VSR
475*
476 IF( ilvsl ) THEN
477 CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
478 $ rwork( iright ), n, vsl, ldvsl, iinfo )
479 IF( iinfo.NE.0 ) THEN
480 info = n + 7
481 GO TO 10
482 END IF
483 END IF
484 IF( ilvsr ) THEN
485 CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
486 $ rwork( iright ), n, vsr, ldvsr, iinfo )
487 IF( iinfo.NE.0 ) THEN
488 info = n + 8
489 GO TO 10
490 END IF
491 END IF
492*
493* Undo scaling
494*
495 IF( ilascl ) THEN
496 CALL zlascl( 'U', -1, -1, anrmto, anrm, n, n, a, lda, iinfo )
497 IF( iinfo.NE.0 ) THEN
498 info = n + 9
499 RETURN
500 END IF
501 CALL zlascl( 'G', -1, -1, anrmto, anrm, n, 1, alpha, n, iinfo )
502 IF( iinfo.NE.0 ) THEN
503 info = n + 9
504 RETURN
505 END IF
506 END IF
507*
508 IF( ilbscl ) THEN
509 CALL zlascl( 'U', -1, -1, bnrmto, bnrm, n, n, b, ldb, iinfo )
510 IF( iinfo.NE.0 ) THEN
511 info = n + 9
512 RETURN
513 END IF
514 CALL zlascl( 'G', -1, -1, bnrmto, bnrm, n, 1, beta, n, iinfo )
515 IF( iinfo.NE.0 ) THEN
516 info = n + 9
517 RETURN
518 END IF
519 END IF
520*
521 10 CONTINUE
522 work( 1 ) = lwkopt
523*
524 RETURN
525*
526* End of ZGEGS
527*
alpha
#define alpha
Definition eval.h:35
zggbak
subroutine zggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
ZGGBAK
Definition zggbak.f:148
zggbal
subroutine zggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
ZGGBAL
Definition zggbal.f:177
zhgeqz
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:284
zlaset
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:106
zgghrd
subroutine zgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
ZGGHRD
Definition zgghrd.f:204
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:128
zgeqrf
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition zgeqrf.f:151

◆ zgegv()

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  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.
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  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 280 of file zgegv.f.

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

◆ zgges()

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.

Definition at line 267 of file zgges.f.

270*
271* -- LAPACK driver routine --
272* -- LAPACK is a software package provided by Univ. of Tennessee, --
273* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274*
275* .. Scalar Arguments ..
276 CHARACTER JOBVSL, JOBVSR, SORT
277 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
278* ..
279* .. Array Arguments ..
280 LOGICAL BWORK( * )
281 DOUBLE PRECISION RWORK( * )
282 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
283 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
284 $ WORK( * )
285* ..
286* .. Function Arguments ..
287 LOGICAL SELCTG
288 EXTERNAL selctg
289* ..
290*
291* =====================================================================
292*
293* .. Parameters ..
294 DOUBLE PRECISION ZERO, ONE
295 parameter( zero = 0.0d0, one = 1.0d0 )
296 COMPLEX*16 CZERO, CONE
297 parameter( czero = ( 0.0d0, 0.0d0 ),
298 $ cone = ( 1.0d0, 0.0d0 ) )
299* ..
300* .. Local Scalars ..
301 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
302 $ LQUERY, WANTST
303 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
304 $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKMIN,
305 $ LWKOPT
306 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
307 $ PVSR, SMLNUM
308* ..
309* .. Local Arrays ..
310 INTEGER IDUM( 1 )
311 DOUBLE PRECISION DIF( 2 )
312* ..
313* .. External Subroutines ..
314 EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghrd,
315 $ zhgeqz, zlacpy, zlascl, zlaset, ztgsen, zungqr,
316 $ zunmqr
317* ..
318* .. External Functions ..
319 LOGICAL LSAME
320 INTEGER ILAENV
321 DOUBLE PRECISION DLAMCH, ZLANGE
322 EXTERNAL lsame, ilaenv, 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 END IF
376*
377* Compute workspace
378* (Note: Comments in the code beginning "Workspace:" describe the
379* minimal amount of workspace needed at that point in the code,
380* as well as the preferred amount for good performance.
381* NB refers to the optimal block size for the immediately
382* following subroutine, as returned by ILAENV.)
383*
384 IF( info.EQ.0 ) THEN
385 lwkmin = max( 1, 2*n )
386 lwkopt = max( 1, n + n*ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
387 lwkopt = max( lwkopt, n +
388 $ n*ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, -1 ) )
389 IF( ilvsl ) THEN
390 lwkopt = max( lwkopt, n +
391 $ n*ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, -1 ) )
392 END IF
393 work( 1 ) = lwkopt
394*
395 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
396 $ info = -18
397 END IF
398*
399 IF( info.NE.0 ) THEN
400 CALL xerbla( 'ZGGES ', -info )
401 RETURN
402 ELSE IF( lquery ) THEN
403 RETURN
404 END IF
405*
406* Quick return if possible
407*
408 IF( n.EQ.0 ) THEN
409 sdim = 0
410 RETURN
411 END IF
412*
413* Get machine constants
414*
415 eps = dlamch( 'P' )
416 smlnum = dlamch( 'S' )
417 bignum = one / smlnum
418 CALL dlabad( smlnum, bignum )
419 smlnum = sqrt( smlnum ) / eps
420 bignum = one / smlnum
421*
422* Scale A if max element outside range [SMLNUM,BIGNUM]
423*
424 anrm = zlange( 'M', n, n, a, lda, rwork )
425 ilascl = .false.
426 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
427 anrmto = smlnum
428 ilascl = .true.
429 ELSE IF( anrm.GT.bignum ) THEN
430 anrmto = bignum
431 ilascl = .true.
432 END IF
433*
434 IF( ilascl )
435 $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
436*
437* Scale B if max element outside range [SMLNUM,BIGNUM]
438*
439 bnrm = zlange( 'M', n, n, b, ldb, rwork )
440 ilbscl = .false.
441 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
442 bnrmto = smlnum
443 ilbscl = .true.
444 ELSE IF( bnrm.GT.bignum ) THEN
445 bnrmto = bignum
446 ilbscl = .true.
447 END IF
448*
449 IF( ilbscl )
450 $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
451*
452* Permute the matrix to make it more nearly triangular
453* (Real Workspace: need 6*N)
454*
455 ileft = 1
456 iright = n + 1
457 irwrk = iright + n
458 CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
459 $ rwork( iright ), rwork( irwrk ), ierr )
460*
461* Reduce B to triangular form (QR decomposition of B)
462* (Complex Workspace: need N, prefer N*NB)
463*
464 irows = ihi + 1 - ilo
465 icols = n + 1 - ilo
466 itau = 1
467 iwrk = itau + irows
468 CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
469 $ work( iwrk ), lwork+1-iwrk, ierr )
470*
471* Apply the orthogonal transformation to matrix A
472* (Complex Workspace: need N, prefer N*NB)
473*
474 CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
475 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
476 $ lwork+1-iwrk, ierr )
477*
478* Initialize VSL
479* (Complex Workspace: need N, prefer N*NB)
480*
481 IF( ilvsl ) THEN
482 CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
483 IF( irows.GT.1 ) THEN
484 CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
485 $ vsl( ilo+1, ilo ), ldvsl )
486 END IF
487 CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
488 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
489 END IF
490*
491* Initialize VSR
492*
493 IF( ilvsr )
494 $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
495*
496* Reduce to generalized Hessenberg form
497* (Workspace: none needed)
498*
499 CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
500 $ ldvsl, vsr, ldvsr, ierr )
501*
502 sdim = 0
503*
504* Perform QZ algorithm, computing Schur vectors if desired
505* (Complex Workspace: need N)
506* (Real Workspace: need N)
507*
508 iwrk = itau
509 CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
510 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
511 $ lwork+1-iwrk, rwork( irwrk ), ierr )
512 IF( ierr.NE.0 ) THEN
513 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
514 info = ierr
515 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
516 info = ierr - n
517 ELSE
518 info = n + 1
519 END IF
520 GO TO 30
521 END IF
522*
523* Sort eigenvalues ALPHA/BETA if desired
524* (Workspace: none needed)
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* (Workspace: none needed)
551*
552 IF( ilvsl )
553 $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
554 $ rwork( iright ), n, vsl, ldvsl, ierr )
555 IF( ilvsr )
556 $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
557 $ rwork( iright ), n, vsr, ldvsr, ierr )
558*
559* Undo scaling
560*
561 IF( ilascl ) THEN
562 CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
563 CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
564 END IF
565*
566 IF( ilbscl ) THEN
567 CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
568 CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
569 END IF
570*
571 IF( wantst ) THEN
572*
573* Check if reordering is correct
574*
575 lastsl = .true.
576 sdim = 0
577 DO 20 i = 1, n
578 cursl = selctg( alpha( i ), beta( i ) )
579 IF( cursl )
580 $ sdim = sdim + 1
581 IF( cursl .AND. .NOT.lastsl )
582 $ info = n + 2
583 lastsl = cursl
584 20 CONTINUE
585*
586 END IF
587*
588 30 CONTINUE
589*
590 work( 1 ) = lwkopt
591*
592 RETURN
593*
594* End of ZGGES
595*
ztgsen
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:433

◆ zgges3()

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 ZLAQZ0
!>                =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.

Definition at line 266 of file zgges3.f.

269*
270* -- LAPACK driver routine --
271* -- LAPACK is a software package provided by Univ. of Tennessee, --
272* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
273*
274* .. Scalar Arguments ..
275 CHARACTER JOBVSL, JOBVSR, SORT
276 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
277* ..
278* .. Array Arguments ..
279 LOGICAL BWORK( * )
280 DOUBLE PRECISION RWORK( * )
281 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
282 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
283 $ WORK( * )
284* ..
285* .. Function Arguments ..
286 LOGICAL SELCTG
287 EXTERNAL selctg
288* ..
289*
290* =====================================================================
291*
292* .. Parameters ..
293 DOUBLE PRECISION ZERO, ONE
294 parameter( zero = 0.0d0, one = 1.0d0 )
295 COMPLEX*16 CZERO, CONE
296 parameter( czero = ( 0.0d0, 0.0d0 ),
297 $ cone = ( 1.0d0, 0.0d0 ) )
298* ..
299* .. Local Scalars ..
300 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
301 $ LQUERY, WANTST
302 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
303 $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT
304 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
305 $ PVSR, SMLNUM
306* ..
307* .. Local Arrays ..
308 INTEGER IDUM( 1 )
309 DOUBLE PRECISION DIF( 2 )
310* ..
311* .. External Subroutines ..
312 EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghd3,
313 $ zlaqz0, zlacpy, zlascl, zlaset, ztgsen, zungqr,
314 $ zunmqr
315* ..
316* .. External Functions ..
317 LOGICAL LSAME
318 DOUBLE PRECISION DLAMCH, ZLANGE
319 EXTERNAL lsame, dlamch, zlange
320* ..
321* .. Intrinsic Functions ..
322 INTRINSIC max, sqrt
323* ..
324* .. Executable Statements ..
325*
326* Decode the input arguments
327*
328 IF( lsame( jobvsl, 'N' ) ) THEN
329 ijobvl = 1
330 ilvsl = .false.
331 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
332 ijobvl = 2
333 ilvsl = .true.
334 ELSE
335 ijobvl = -1
336 ilvsl = .false.
337 END IF
338*
339 IF( lsame( jobvsr, 'N' ) ) THEN
340 ijobvr = 1
341 ilvsr = .false.
342 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
343 ijobvr = 2
344 ilvsr = .true.
345 ELSE
346 ijobvr = -1
347 ilvsr = .false.
348 END IF
349*
350 wantst = lsame( sort, 'S' )
351*
352* Test the input arguments
353*
354 info = 0
355 lquery = ( lwork.EQ.-1 )
356 IF( ijobvl.LE.0 ) THEN
357 info = -1
358 ELSE IF( ijobvr.LE.0 ) THEN
359 info = -2
360 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
361 info = -3
362 ELSE IF( n.LT.0 ) THEN
363 info = -5
364 ELSE IF( lda.LT.max( 1, n ) ) THEN
365 info = -7
366 ELSE IF( ldb.LT.max( 1, n ) ) THEN
367 info = -9
368 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
369 info = -14
370 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
371 info = -16
372 ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
373 info = -18
374 END IF
375*
376* Compute workspace
377*
378 IF( info.EQ.0 ) THEN
379 CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
380 lwkopt = max( 1, n + int( work( 1 ) ) )
381 CALL zunmqr( 'L', 'C', n, n, n, b, ldb, work, a, lda, work,
382 $ -1, ierr )
383 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
384 IF( ilvsl ) THEN
385 CALL zungqr( n, n, n, vsl, ldvsl, work, work, -1, ierr )
386 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
387 END IF
388 CALL zgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
389 $ ldvsl, vsr, ldvsr, work, -1, ierr )
390 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
391 CALL zlaqz0( 'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
392 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work, -1,
393 $ rwork, 0, ierr )
394 lwkopt = max( lwkopt, int( work( 1 ) ) )
395 IF( wantst ) THEN
396 CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
397 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, sdim,
398 $ pvsl, pvsr, dif, work, -1, idum, 1, ierr )
399 lwkopt = max( lwkopt, int( work( 1 ) ) )
400 END IF
401 work( 1 ) = dcmplx( lwkopt )
402 END IF
403*
404 IF( info.NE.0 ) THEN
405 CALL xerbla( 'ZGGES3 ', -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 ) THEN
414 sdim = 0
415 RETURN
416 END IF
417*
418* Get machine constants
419*
420 eps = dlamch( 'P' )
421 smlnum = dlamch( 'S' )
422 bignum = one / smlnum
423 CALL dlabad( smlnum, bignum )
424 smlnum = sqrt( smlnum ) / eps
425 bignum = one / smlnum
426*
427* Scale A if max element outside range [SMLNUM,BIGNUM]
428*
429 anrm = zlange( 'M', n, n, a, lda, rwork )
430 ilascl = .false.
431 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
432 anrmto = smlnum
433 ilascl = .true.
434 ELSE IF( anrm.GT.bignum ) THEN
435 anrmto = bignum
436 ilascl = .true.
437 END IF
438*
439 IF( ilascl )
440 $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
441*
442* Scale B if max element outside range [SMLNUM,BIGNUM]
443*
444 bnrm = zlange( 'M', n, n, b, ldb, rwork )
445 ilbscl = .false.
446 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
447 bnrmto = smlnum
448 ilbscl = .true.
449 ELSE IF( bnrm.GT.bignum ) THEN
450 bnrmto = bignum
451 ilbscl = .true.
452 END IF
453*
454 IF( ilbscl )
455 $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
456*
457* Permute the matrix to make it more nearly triangular
458*
459 ileft = 1
460 iright = n + 1
461 irwrk = iright + n
462 CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
463 $ rwork( iright ), rwork( irwrk ), ierr )
464*
465* Reduce B to triangular form (QR decomposition of B)
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*
476 CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
477 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
478 $ lwork+1-iwrk, ierr )
479*
480* Initialize VSL
481*
482 IF( ilvsl ) THEN
483 CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
484 IF( irows.GT.1 ) THEN
485 CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
486 $ vsl( ilo+1, ilo ), ldvsl )
487 END IF
488 CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
489 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
490 END IF
491*
492* Initialize VSR
493*
494 IF( ilvsr )
495 $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
496*
497* Reduce to generalized Hessenberg form
498*
499 CALL zgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
500 $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk, ierr )
501*
502 sdim = 0
503*
504* Perform QZ algorithm, computing Schur vectors if desired
505*
506 iwrk = itau
507 CALL zlaqz0( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
508 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
509 $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
510 IF( ierr.NE.0 ) THEN
511 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
512 info = ierr
513 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
514 info = ierr - n
515 ELSE
516 info = n + 1
517 END IF
518 GO TO 30
519 END IF
520*
521* Sort eigenvalues ALPHA/BETA if desired
522*
523 IF( wantst ) THEN
524*
525* Undo scaling on eigenvalues before selecting
526*
527 IF( ilascl )
528 $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
529 IF( ilbscl )
530 $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, 1, beta, n, ierr )
531*
532* Select eigenvalues
533*
534 DO 10 i = 1, n
535 bwork( i ) = selctg( alpha( i ), beta( i ) )
536 10 CONTINUE
537*
538 CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alpha,
539 $ beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl, pvsr,
540 $ dif, work( iwrk ), lwork-iwrk+1, idum, 1, ierr )
541 IF( ierr.EQ.1 )
542 $ info = n + 3
543*
544 END IF
545*
546* Apply back-permutation to VSL and VSR
547*
548 IF( ilvsl )
549 $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
550 $ rwork( iright ), n, vsl, ldvsl, ierr )
551 IF( ilvsr )
552 $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
553 $ rwork( iright ), n, vsr, ldvsr, ierr )
554*
555* Undo scaling
556*
557 IF( ilascl ) THEN
558 CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
559 CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
560 END IF
561*
562 IF( ilbscl ) THEN
563 CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
564 CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
565 END IF
566*
567 IF( wantst ) THEN
568*
569* Check if reordering is correct
570*
571 lastsl = .true.
572 sdim = 0
573 DO 20 i = 1, n
574 cursl = selctg( alpha( i ), beta( i ) )
575 IF( cursl )
576 $ sdim = sdim + 1
577 IF( cursl .AND. .NOT.lastsl )
578 $ info = n + 2
579 lastsl = cursl
580 20 CONTINUE
581*
582 END IF
583*
584 30 CONTINUE
585*
586 work( 1 ) = dcmplx( lwkopt )
587*
588 RETURN
589*
590* End of ZGGES3
591*
zlaqz0
recursive subroutine zlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, rec, info)
ZLAQZ0
Definition zlaqz0.f:284
zgghd3
subroutine zgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
ZGGHD3
Definition zgghd3.f:227

◆ zggesx()

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 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.
!>          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.

Definition at line 326 of file zggesx.f.

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

◆ zggev()

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 ZHGEQZ,
!>                =N+2: error return from ZTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 215 of file zggev.f.

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

◆ zggev3()

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 ZHGEQZ,
!>                =N+2: error return from ZTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 214 of file zggev3.f.

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

◆ zggevx()

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  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  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.
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 370 of file zggevx.f.

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