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

Functions

subroutine cgees (jobvs, sort, select, n, a, lda, sdim, w, vs, ldvs, work, lwork, rwork, bwork, info)
  CGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine cgeesx (jobvs, sort, select, sense, n, a, lda, sdim, w, vs, ldvs, rconde, rcondv, work, lwork, rwork, bwork, info)
  CGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine cgeev (jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
  CGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine cgeevx (balanc, jobvl, jobvr, sense, n, a, lda, w, vl, ldvl, vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv, work, lwork, rwork, info)
  CGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine cgges (jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alpha, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, rwork, bwork, info)
  CGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine cgges3 (jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alpha, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, rwork, bwork, info)
  CGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)
subroutine cggesx (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)
  CGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine cggev (jobvl, jobvr, n, a, lda, b, ldb, alpha, beta, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
  CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine cggev3 (jobvl, jobvr, n, a, lda, b, ldb, alpha, beta, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
  CGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)
subroutine cggevx (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)
  CGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine cgegs (jobvsl, jobvsr, n, a, lda, b, ldb, alpha, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, rwork, info)
  CGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine cgegv (jobvl, jobvr, n, a, lda, b, ldb, alpha, beta, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
  CGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

Detailed Description

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

Function Documentation

◆ cgees()

subroutine cgees ( character jobvs,
character sort,
external select,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer sdim,
complex, dimension( * ) w,
complex, dimension( ldvs, * ) vs,
integer ldvs,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> CGEES 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 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 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 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 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 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 REAL 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 cgees.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 REAL RWORK( * )
209 COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
210* ..
211* .. Function Arguments ..
212 LOGICAL SELECT
213 EXTERNAL SELECT
214* ..
215*
216* =====================================================================
217*
218* .. Parameters ..
219 REAL ZERO, ONE
220 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
227* ..
228* .. Local Arrays ..
229 REAL DUM( 1 )
230* ..
231* .. External Subroutines ..
232 EXTERNAL ccopy, cgebak, cgebal, cgehrd, chseqr, clacpy,
233 $ clascl, ctrsen, cunghr, slabad, xerbla
234* ..
235* .. External Functions ..
236 LOGICAL LSAME
237 INTEGER ILAENV
238 REAL CLANGE, SLAMCH
239 EXTERNAL lsame, ilaenv, clange, slamch
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 CHSEQR, 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, 'CGEHRD', ' ', n, 1, n, 0 )
281 minwrk = 2*n
282*
283 CALL chseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
284 $ work, -1, ieval )
285 hswork = real( work( 1 ) )
286*
287 IF( .NOT.wantvs ) THEN
288 maxwrk = max( maxwrk, hswork )
289 ELSE
290 maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'CUNGHR',
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( 'CGEES ', -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 = slamch( 'P' )
319 smlnum = slamch( 'S' )
320 bignum = one / smlnum
321 CALL slabad( 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 = clange( '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 clascl( '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 cgebal( '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 cgehrd( 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 clacpy( '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 cunghr( 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 chseqr( '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 clascl( '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 ctrsen( '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 cgebak( '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 clascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
413 CALL ccopy( n, a, lda+1, w, 1 )
414 END IF
415*
416 work( 1 ) = maxwrk
417 RETURN
418*
419* End of CGEES
420*
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.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
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
cgebal
subroutine cgebal(job, n, a, lda, ilo, ihi, scale, info)
CGEBAL
Definition cgebal.f:161
cgehrd
subroutine cgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
CGEHRD
Definition cgehrd.f:167
cgebak
subroutine cgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
CGEBAK
Definition cgebak.f:131
clascl
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
cunghr
subroutine cunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
CUNGHR
Definition cunghr.f:126
ctrsen
subroutine ctrsen(job, compq, select, n, t, ldt, q, ldq, w, m, s, sep, work, lwork, info)
CTRSEN
Definition ctrsen.f:264
chseqr
subroutine chseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
CHSEQR
Definition chseqr.f:299
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
max
#define max(a, b)
Definition macros.h:21

◆ cgeesx()

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

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

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

Purpose:
!>
!> CGEESX 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 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 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 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 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 REAL
!>          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 REAL
!>          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 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 REAL 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 cgeesx.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 REAL RCONDE, RCONDV
248* ..
249* .. Array Arguments ..
250 LOGICAL BWORK( * )
251 REAL RWORK( * )
252 COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
253* ..
254* .. Function Arguments ..
255 LOGICAL SELECT
256 EXTERNAL SELECT
257* ..
258*
259* =====================================================================
260*
261* .. Parameters ..
262 REAL ZERO, ONE
263 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, CSCALE, EPS, SMLNUM
271* ..
272* .. Local Arrays ..
273 REAL DUM( 1 )
274* ..
275* .. External Subroutines ..
276 EXTERNAL ccopy, cgebak, cgebal, cgehrd, chseqr, clacpy,
277 $ clascl, ctrsen, cunghr, slabad, slascl, xerbla
278* ..
279* .. External Functions ..
280 LOGICAL LSAME
281 INTEGER ILAENV
282 REAL CLANGE, SLAMCH
283 EXTERNAL lsame, ilaenv, clange, slamch
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 CHSEQR, 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 CTRSEN 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, 'CGEHRD', ' ', n, 1, n, 0 )
336 minwrk = 2*n
337*
338 CALL chseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
339 $ work, -1, ieval )
340 hswork = real( work( 1 ) )
341*
342 IF( .NOT.wantvs ) THEN
343 maxwrk = max( maxwrk, hswork )
344 ELSE
345 maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'CUNGHR',
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( 'CGEESX', -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 = slamch( 'P' )
377 smlnum = slamch( 'S' )
378 bignum = one / smlnum
379 CALL slabad( 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 = clange( '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 clascl( '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 cgebal( '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 cgehrd( 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 clacpy( '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 cunghr( 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 chseqr( '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 clascl( '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 ctrsen( 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 cgebak( '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 clascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
483 CALL ccopy( n, a, lda+1, w, 1 )
484 IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
485 dum( 1 ) = rcondv
486 CALL slascl( '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 CGEESX
495*
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:143

◆ cgeev()

subroutine cgeev ( character jobvl,
character jobvr,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( * ) w,
complex, dimension( ldvl, * ) vl,
integer ldvl,
complex, dimension( ldvr, * ) vr,
integer ldvr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CGEEV 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 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 array, dimension (N)
!>          W contains the computed eigenvalues.
!>
[out]VL
!>          VL is COMPLEX 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 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 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 REAL 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 cgeev.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 REAL RWORK( * )
192 COMPLEX A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
193 $ W( * ), WORK( * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO, ONE
200 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
208 COMPLEX TMP
209* ..
210* .. Local Arrays ..
211 LOGICAL SELECT( 1 )
212 REAL DUM( 1 )
213* ..
214* .. External Subroutines ..
215 EXTERNAL slabad, xerbla, csscal, cgebak, cgebal, cgehrd,
216 $ chseqr, clacpy, clascl, cscal, ctrevc3, cunghr
217* ..
218* .. External Functions ..
219 LOGICAL LSAME
220 INTEGER ISAMAX, ILAENV
221 REAL SLAMCH, SCNRM2, CLANGE
222 EXTERNAL lsame, isamax, ilaenv, slamch, scnrm2, clange
223* ..
224* .. Intrinsic Functions ..
225 INTRINSIC real, cmplx, 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 CHSEQR, 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, 'CGEHRD', ' ', n, 1, n, 0 )
266 minwrk = 2*n
267 IF( wantvl ) THEN
268 maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'CUNGHR',
269 $ ' ', n, 1, n, -1 ) )
270 CALL ctrevc3( '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 chseqr( '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, 'CUNGHR',
279 $ ' ', n, 1, n, -1 ) )
280 CALL ctrevc3( '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 chseqr( 'S', 'V', n, 1, n, a, lda, w, vr, ldvr,
286 $ work, -1, info )
287 ELSE
288 CALL chseqr( '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( 'CGEEV ', -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 = slamch( 'P' )
316 smlnum = slamch( 'S' )
317 bignum = one / smlnum
318 CALL slabad( 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 = clange( '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 clascl( '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 cgebal( '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 cgehrd( 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 clacpy( '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 cunghr( 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 chseqr( '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 clacpy( '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 clacpy( '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 cunghr( 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 chseqr( '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 chseqr( '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 CHSEQR, 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 ctrevc3( 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 cgebak( '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 / scnrm2( n, vl( 1, i ), 1 )
448 CALL csscal( n, scl, vl( 1, i ), 1 )
449 DO 10 k = 1, n
450 rwork( irwork+k-1 ) = real( vl( k, i ) )**2 +
451 $ aimag( vl( k, i ) )**2
452 10 CONTINUE
453 k = isamax( n, rwork( irwork ), 1 )
454 tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
455 CALL cscal( n, tmp, vl( 1, i ), 1 )
456 vl( k, i ) = cmplx( real( 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 cgebak( '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 / scnrm2( n, vr( 1, i ), 1 )
473 CALL csscal( n, scl, vr( 1, i ), 1 )
474 DO 30 k = 1, n
475 rwork( irwork+k-1 ) = real( vr( k, i ) )**2 +
476 $ aimag( vr( k, i ) )**2
477 30 CONTINUE
478 k = isamax( n, rwork( irwork ), 1 )
479 tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
480 CALL cscal( n, tmp, vr( 1, i ), 1 )
481 vr( k, i ) = cmplx( real( 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 clascl( '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 clascl( '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 CGEEV
500*
cmplx
float cmplx[2]
Definition pblas.h:136
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
ctrevc3
subroutine ctrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, rwork, lrwork, info)
CTREVC3
Definition ctrevc3.f:244
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
scnrm2
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90

◆ cgeevx()

subroutine cgeevx ( character balanc,
character jobvl,
character jobvr,
character sense,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( * ) w,
complex, dimension( ldvl, * ) vl,
integer ldvl,
complex, dimension( ldvr, * ) vr,
integer ldvr,
integer ilo,
integer ihi,
real, dimension( * ) scale,
real abnrm,
real, dimension( * ) rconde,
real, dimension( * ) rcondv,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CGEEVX 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 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 array, dimension (N)
!>          W contains the computed eigenvalues.
!>
[out]VL
!>          VL is COMPLEX 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 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 REAL 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 REAL
!>          The one-norm of the balanced matrix (the maximum
!>          of the sum of absolute values of elements of any column).
!>
[out]RCONDE
!>          RCONDE is REAL array, dimension (N)
!>          RCONDE(j) is the reciprocal condition number of the j-th
!>          eigenvalue.
!>
[out]RCONDV
!>          RCONDV is REAL array, dimension (N)
!>          RCONDV(j) is the reciprocal condition number of the j-th
!>          right eigenvector.
!>
[out]WORK
!>          WORK is COMPLEX 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 REAL 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 cgeevx.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 REAL ABNRM
298* ..
299* .. Array Arguments ..
300 REAL RCONDE( * ), RCONDV( * ), RWORK( * ),
301 $ SCALE( * )
302 COMPLEX A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
303 $ W( * ), WORK( * )
304* ..
305*
306* =====================================================================
307*
308* .. Parameters ..
309 REAL ZERO, ONE
310 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
319 COMPLEX TMP
320* ..
321* .. Local Arrays ..
322 LOGICAL SELECT( 1 )
323 REAL DUM( 1 )
324* ..
325* .. External Subroutines ..
326 EXTERNAL slabad, slascl, xerbla, csscal, cgebak, cgebal,
327 $ cgehrd, chseqr, clacpy, clascl, cscal, ctrevc3,
328 $ ctrsna, cunghr
329* ..
330* .. External Functions ..
331 LOGICAL LSAME
332 INTEGER ISAMAX, ILAENV
333 REAL SLAMCH, SCNRM2, CLANGE
334 EXTERNAL lsame, isamax, ilaenv, slamch, scnrm2, clange
335* ..
336* .. Intrinsic Functions ..
337 INTRINSIC real, cmplx, 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 CHSEQR, 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, 'CGEHRD', ' ', n, 1, n, 0 )
389*
390 IF( wantvl ) THEN
391 CALL ctrevc3( '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 chseqr( 'S', 'V', n, 1, n, a, lda, w, vl, ldvl,
397 $ work, -1, info )
398 ELSE IF( wantvr ) THEN
399 CALL ctrevc3( '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 chseqr( 'S', 'V', n, 1, n, a, lda, w, vr, ldvr,
405 $ work, -1, info )
406 ELSE
407 IF( wntsnn ) THEN
408 CALL chseqr( 'E', 'N', n, 1, n, a, lda, w, vr, ldvr,
409 $ work, -1, info )
410 ELSE
411 CALL chseqr( '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, 'CUNGHR',
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( 'CGEEVX', -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 = slamch( 'P' )
459 smlnum = slamch( 'S' )
460 bignum = one / smlnum
461 CALL slabad( 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 = clange( '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 clascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
479*
480* Balance the matrix and compute ABNRM
481*
482 CALL cgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
483 abnrm = clange( '1', n, n, a, lda, dum )
484 IF( scalea ) THEN
485 dum( 1 ) = abnrm
486 CALL slascl( '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 cgehrd( 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 clacpy( '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 cunghr( 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 chseqr( '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 clacpy( '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 clacpy( '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 cunghr( 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 chseqr( '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 chseqr( 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 CHSEQR, 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 ctrevc3( 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 ctrsna( 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 cgebak( 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 / scnrm2( n, vl( 1, i ), 1 )
610 CALL csscal( n, scl, vl( 1, i ), 1 )
611 DO 10 k = 1, n
612 rwork( k ) = real( vl( k, i ) )**2 +
613 $ aimag( vl( k, i ) )**2
614 10 CONTINUE
615 k = isamax( n, rwork, 1 )
616 tmp = conjg( vl( k, i ) ) / sqrt( rwork( k ) )
617 CALL cscal( n, tmp, vl( 1, i ), 1 )
618 vl( k, i ) = cmplx( real( 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 cgebak( 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 / scnrm2( n, vr( 1, i ), 1 )
633 CALL csscal( n, scl, vr( 1, i ), 1 )
634 DO 30 k = 1, n
635 rwork( k ) = real( vr( k, i ) )**2 +
636 $ aimag( vr( k, i ) )**2
637 30 CONTINUE
638 k = isamax( n, rwork, 1 )
639 tmp = conjg( vr( k, i ) ) / sqrt( rwork( k ) )
640 CALL cscal( n, tmp, vr( 1, i ), 1 )
641 vr( k, i ) = cmplx( real( 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 clascl( '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 slascl( 'G', 0, 0, cscale, anrm, n, 1, rcondv, n,
654 $ ierr )
655 ELSE
656 CALL clascl( '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 CGEEVX
664*
ctrsna
subroutine ctrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, rwork, info)
CTRSNA
Definition ctrsna.f:249

◆ cgegs()

subroutine cgegs ( character jobvsl,
character jobvsr,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvsl, * ) vsl,
integer ldvsl,
complex, dimension( ldvsr, * ) vsr,
integer ldvsr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine CGGES.
!>
!> CGEGS 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
!> CGEGV should be used instead.  See CGEGV 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 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 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 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 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 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 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 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 CGEQRF, CUNMQR, and CUNGQR.)  Then compute:
!>          NB  -- MAX of the blocksizes for CGEQRF, CUNMQR, 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 REAL 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 CGGBAL
!>                =N+2: error return from CGEQRF
!>                =N+3: error return from CUNMQR
!>                =N+4: error return from CUNGQR
!>                =N+5: error return from CGGHRD
!>                =N+6: error return from CHGEQZ (other than failed
!>                                               iteration)
!>                =N+7: error return from CGGBAK (computing VSL)
!>                =N+8: error return from CGGBAK (computing VSR)
!>                =N+9: error return from CLASCL (various places)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 222 of file cgegs.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 REAL RWORK( * )
236 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
237 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
238 $ WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 REAL ZERO, ONE
245 parameter( zero = 0.0e0, one = 1.0e0 )
246 COMPLEX CZERO, CONE
247 parameter( czero = ( 0.0e0, 0.0e0 ),
248 $ cone = ( 1.0e0, 0.0e0 ) )
249* ..
250* .. Local Scalars ..
251 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
252 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT,
253 $ ILO, IRIGHT, IROWS, IRWORK, ITAU, IWORK,
254 $ LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
255 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
256 $ SAFMIN, SMLNUM
257* ..
258* .. External Subroutines ..
259 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
260 $ clascl, claset, cungqr, cunmqr, xerbla
261* ..
262* .. External Functions ..
263 LOGICAL LSAME
264 INTEGER ILAENV
265 REAL CLANGE, SLAMCH
266 EXTERNAL ilaenv, lsame, clange, slamch
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, 'CGEQRF', ' ', n, n, -1, -1 )
324 nb2 = ilaenv( 1, 'CUNMQR', ' ', n, n, n, -1 )
325 nb3 = ilaenv( 1, 'CUNGQR', ' ', 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( 'CGEGS ', -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 = slamch( 'E' )*slamch( 'B' )
346 safmin = slamch( 'S' )
347 smlnum = n*safmin / eps
348 bignum = one / smlnum
349*
350* Scale A if max element outside range [SMLNUM,BIGNUM]
351*
352 anrm = clange( '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 clascl( '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 = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vsl, ldvsl )
430 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
431 $ vsl( ilo+1, ilo ), ldvsl )
432 CALL cungqr( 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 claset( 'Full', n, n, czero, cone, vsr, ldvsr )
445*
446* Reduce to generalized Hessenberg form
447*
448 CALL cgghrd( 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 chgeqz( '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 cggbak( '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 cggbak( '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 clascl( '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 clascl( '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 clascl( '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 clascl( '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 CGEGS
527*
alpha
#define alpha
Definition eval.h:35
cggbak
subroutine cggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
CGGBAK
Definition cggbak.f:148
cggbal
subroutine cggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
CGGBAL
Definition cggbal.f:177
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
chgeqz
subroutine chgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, info)
CHGEQZ
Definition chgeqz.f:284
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
cunmqr
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
cgghrd
subroutine cgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
CGGHRD
Definition cgghrd.f:204
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128

◆ cgegv()

subroutine cgegv ( character jobvl,
character jobvr,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvl, * ) vl,
integer ldvl,
complex, dimension( ldvr, * ) vr,
integer ldvr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine CGGEV.
!>
!> CGEGV 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 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 CGGHRD and CHGEQZ
!>          for details.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is COMPLEX 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 CGGHRD and CHGEQZ for
!>          details.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHA
!>          ALPHA is COMPLEX array, dimension (N)
!>          The complex scalars alpha that define the eigenvalues of
!>          GNEP.
!>
[out]BETA
!>          BETA is COMPLEX 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 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 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 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 CGEQRF, CUNMQR, and CUNGQR.)  Then compute:
!>          NB  -- MAX of the blocksizes for CGEQRF, CUNMQR, and CUNGQR;
!>          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 REAL 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 CGGBAL
!>                =N+2: error return from CGEQRF
!>                =N+3: error return from CUNMQR
!>                =N+4: error return from CUNGQR
!>                =N+5: error return from CGGHRD
!>                =N+6: error return from CHGEQZ (other than failed
!>                                               iteration)
!>                =N+7: error return from CTGEVC
!>                =N+8: error return from CGGBAK (computing VL)
!>                =N+9: error return from CGGBAK (computing VR)
!>                =N+10: error return from CLASCL (various calls)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Balancing
!>  ---------
!>
!>  This driver calls CGGBAL 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, CGGBAK 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 cgegv.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 REAL RWORK( * )
293 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
294 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
295 $ WORK( * )
296* ..
297*
298* =====================================================================
299*
300* .. Parameters ..
301 REAL ZERO, ONE
302 parameter( zero = 0.0e0, one = 1.0e0 )
303 COMPLEX CZERO, CONE
304 parameter( czero = ( 0.0e0, 0.0e0 ),
305 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
314 $ BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
315 $ SALFAR, SBETA, SCALE, TEMP
316 COMPLEX X
317* ..
318* .. Local Arrays ..
319 LOGICAL LDUMMA( 1 )
320* ..
321* .. External Subroutines ..
322 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
323 $ clascl, claset, ctgevc, cungqr, cunmqr, xerbla
324* ..
325* .. External Functions ..
326 LOGICAL LSAME
327 INTEGER ILAENV
328 REAL CLANGE, SLAMCH
329 EXTERNAL ilaenv, lsame, clange, slamch
330* ..
331* .. Intrinsic Functions ..
332 INTRINSIC abs, aimag, cmplx, int, max, real
333* ..
334* .. Statement Functions ..
335 REAL ABS1
336* ..
337* .. Statement Function definitions ..
338 abs1( x ) = abs( real( x ) ) + abs( aimag( 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, 'CGEQRF', ' ', n, n, -1, -1 )
394 nb2 = ilaenv( 1, 'CUNMQR', ' ', n, n, n, -1 )
395 nb3 = ilaenv( 1, 'CUNGQR', ' ', 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( 'CGEGV ', -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 = slamch( 'E' )*slamch( 'B' )
416 safmin = slamch( 'S' )
417 safmin = safmin + safmin
418 safmax = one / safmin
419*
420* Scale A
421*
422 anrm = clange( '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 clascl( '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 = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vl, ldvl )
504 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
505 $ vl( ilo+1, ilo ), ldvl )
506 CALL cungqr( 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 claset( '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 cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
527 $ ldvl, vr, ldvr, iinfo )
528 ELSE
529 CALL cgghrd( '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 chgeqz( 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 ctgevc( 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 cggbak( '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 cggbak( '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( real( alpha( jc ) ) )
640 absai = abs( aimag( alpha( jc ) ) )
641 absb = abs( real( beta( jc ) ) )
642 salfar = anrm*real( alpha( jc ) )
643 salfai = anrm*aimag( alpha( jc ) )
644 sbeta = bnrm*real( 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*real( alpha( jc ) ) )*anrm
689 salfai = ( scale*aimag( alpha( jc ) ) )*anrm
690 sbeta = ( scale*beta( jc ) )*bnrm
691 END IF
692 alpha( jc ) = cmplx( salfar, salfai )
693 beta( jc ) = sbeta
694 70 CONTINUE
695*
696 80 CONTINUE
697 work( 1 ) = lwkopt
698*
699 RETURN
700*
701* End of CGEGV
702*
ctgevc
subroutine ctgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
CTGEVC
Definition ctgevc.f:219
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339

◆ cgges()

subroutine cgges ( character jobvsl,
character jobvsr,
character sort,
external selctg,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
integer sdim,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvsl, * ) vsl,
integer ldvsl,
complex, dimension( ldvsr, * ) vsr,
integer ldvsr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> CGGES 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
!> CGGEV 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 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 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 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 array, dimension (N)
!>
[out]BETA
!>          BETA is COMPLEX 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 CGGES. 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 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 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 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 REAL 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 CHGEQZ
!>                =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 CTGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 267 of file cgges.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 REAL RWORK( * )
282 COMPLEX 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 REAL ZERO, ONE
295 parameter( zero = 0.0e0, one = 1.0e0 )
296 COMPLEX CZERO, CONE
297 parameter( czero = ( 0.0e0, 0.0e0 ),
298 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
307 $ PVSR, SMLNUM
308* ..
309* .. Local Arrays ..
310 INTEGER IDUM( 1 )
311 REAL DIF( 2 )
312* ..
313* .. External Subroutines ..
314 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
315 $ clascl, claset, ctgsen, cungqr, cunmqr, slabad,
316 $ xerbla
317* ..
318* .. External Functions ..
319 LOGICAL LSAME
320 INTEGER ILAENV
321 REAL CLANGE, SLAMCH
322 EXTERNAL lsame, ilaenv, clange, slamch
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, 'CGEQRF', ' ', n, 1, n, 0 ) )
387 lwkopt = max( lwkopt, n +
388 $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, -1 ) )
389 IF( ilvsl ) THEN
390 lwkopt = max( lwkopt, n +
391 $ n*ilaenv( 1, 'CUNGQR', ' ', 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( 'CGGES ', -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 = slamch( 'P' )
416 smlnum = slamch( 'S' )
417 bignum = one / smlnum
418 CALL slabad( 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 = clange( '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 clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
436*
437* Scale B if max element outside range [SMLNUM,BIGNUM]
438*
439 bnrm = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vsl, ldvsl )
483 IF( irows.GT.1 ) THEN
484 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
485 $ vsl( ilo+1, ilo ), ldvsl )
486 END IF
487 CALL cungqr( 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 claset( 'Full', n, n, czero, cone, vsr, ldvsr )
495*
496* Reduce to generalized Hessenberg form
497* (Workspace: none needed)
498*
499 CALL cgghrd( 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 chgeqz( '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 clascl( 'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
532 IF( ilbscl )
533 $ CALL clascl( '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 ctgsen( 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 cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
554 $ rwork( iright ), n, vsl, ldvsl, ierr )
555 IF( ilvsr )
556 $ CALL cggbak( '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 clascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
563 CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
564 END IF
565*
566 IF( ilbscl ) THEN
567 CALL clascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
568 CALL clascl( '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 CGGES
595*
ctgsen
subroutine ctgsen(ijob, wantq, wantz, select, n, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, m, pl, pr, dif, work, lwork, iwork, liwork, info)
CTGSEN
Definition ctgsen.f:433

◆ cgges3()

subroutine cgges3 ( character jobvsl,
character jobvsr,
character sort,
external selctg,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
integer sdim,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvsl, * ) vsl,
integer ldvsl,
complex, dimension( ldvsr, * ) vsr,
integer ldvsr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> CGGES3 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
!> CGGEV 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 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 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 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 array, dimension (N)
!>
[out]BETA
!>          BETA is COMPLEX 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 CGGES3. 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 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 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 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 REAL 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 CLAQZ0
!>                =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 CTGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 266 of file cgges3.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 REAL RWORK( * )
281 COMPLEX 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 REAL ZERO, ONE
294 parameter( zero = 0.0e0, one = 1.0e0 )
295 COMPLEX CZERO, CONE
296 parameter( czero = ( 0.0e0, 0.0e0 ),
297 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
305 $ PVSR, SMLNUM
306* ..
307* .. Local Arrays ..
308 INTEGER IDUM( 1 )
309 REAL DIF( 2 )
310* ..
311* .. External Subroutines ..
312 EXTERNAL cgeqrf, cggbak, cggbal, cgghd3, claqz0, clacpy,
313 $ clascl, claset, ctgsen, cungqr, cunmqr, slabad,
314 $ xerbla
315* ..
316* .. External Functions ..
317 LOGICAL LSAME
318 REAL CLANGE, SLAMCH
319 EXTERNAL lsame, clange, slamch
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 cgeqrf( n, n, b, ldb, work, work, -1, ierr )
380 lwkopt = max( 1, n + int( work( 1 ) ) )
381 CALL cunmqr( '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 cungqr( n, n, n, vsl, ldvsl, work, work, -1,
386 $ ierr )
387 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
388 END IF
389 CALL cgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
390 $ ldvsl, vsr, ldvsr, work, -1, ierr )
391 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
392 CALL claqz0( 'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
393 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work, -1,
394 $ rwork, 0, ierr )
395 lwkopt = max( lwkopt, int( work( 1 ) ) )
396 IF( wantst ) THEN
397 CALL ctgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
398 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, sdim,
399 $ pvsl, pvsr, dif, work, -1, idum, 1, ierr )
400 lwkopt = max( lwkopt, int( work( 1 ) ) )
401 END IF
402 work( 1 ) = cmplx( lwkopt )
403 END IF
404
405*
406 IF( info.NE.0 ) THEN
407 CALL xerbla( 'CGGES3 ', -info )
408 RETURN
409 ELSE IF( lquery ) THEN
410 RETURN
411 END IF
412*
413* Quick return if possible
414*
415 IF( n.EQ.0 ) THEN
416 sdim = 0
417 RETURN
418 END IF
419*
420* Get machine constants
421*
422 eps = slamch( 'P' )
423 smlnum = slamch( 'S' )
424 bignum = one / smlnum
425 CALL slabad( smlnum, bignum )
426 smlnum = sqrt( smlnum ) / eps
427 bignum = one / smlnum
428*
429* Scale A if max element outside range [SMLNUM,BIGNUM]
430*
431 anrm = clange( 'M', n, n, a, lda, rwork )
432 ilascl = .false.
433 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
434 anrmto = smlnum
435 ilascl = .true.
436 ELSE IF( anrm.GT.bignum ) THEN
437 anrmto = bignum
438 ilascl = .true.
439 END IF
440*
441 IF( ilascl )
442 $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
443*
444* Scale B if max element outside range [SMLNUM,BIGNUM]
445*
446 bnrm = clange( 'M', n, n, b, ldb, rwork )
447 ilbscl = .false.
448 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
449 bnrmto = smlnum
450 ilbscl = .true.
451 ELSE IF( bnrm.GT.bignum ) THEN
452 bnrmto = bignum
453 ilbscl = .true.
454 END IF
455*
456 IF( ilbscl )
457 $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
458*
459* Permute the matrix to make it more nearly triangular
460*
461 ileft = 1
462 iright = n + 1
463 irwrk = iright + n
464 CALL cggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
465 $ rwork( iright ), rwork( irwrk ), ierr )
466*
467* Reduce B to triangular form (QR decomposition of B)
468*
469 irows = ihi + 1 - ilo
470 icols = n + 1 - ilo
471 itau = 1
472 iwrk = itau + irows
473 CALL cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
474 $ work( iwrk ), lwork+1-iwrk, ierr )
475*
476* Apply the orthogonal transformation to matrix A
477*
478 CALL cunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
479 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
480 $ lwork+1-iwrk, ierr )
481*
482* Initialize VSL
483*
484 IF( ilvsl ) THEN
485 CALL claset( 'Full', n, n, czero, cone, vsl, ldvsl )
486 IF( irows.GT.1 ) THEN
487 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
488 $ vsl( ilo+1, ilo ), ldvsl )
489 END IF
490 CALL cungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
491 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
492 END IF
493*
494* Initialize VSR
495*
496 IF( ilvsr )
497 $ CALL claset( 'Full', n, n, czero, cone, vsr, ldvsr )
498*
499* Reduce to generalized Hessenberg form
500*
501 CALL cgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
502 $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk, ierr )
503*
504 sdim = 0
505*
506* Perform QZ algorithm, computing Schur vectors if desired
507*
508 iwrk = itau
509 CALL claqz0( '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 ), 0, 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*
525 IF( wantst ) THEN
526*
527* Undo scaling on eigenvalues before selecting
528*
529 IF( ilascl )
530 $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
531 IF( ilbscl )
532 $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, 1, beta, n, ierr )
533*
534* Select eigenvalues
535*
536 DO 10 i = 1, n
537 bwork( i ) = selctg( alpha( i ), beta( i ) )
538 10 CONTINUE
539*
540 CALL ctgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alpha,
541 $ beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl, pvsr,
542 $ dif, work( iwrk ), lwork-iwrk+1, idum, 1, ierr )
543 IF( ierr.EQ.1 )
544 $ info = n + 3
545*
546 END IF
547*
548* Apply back-permutation to VSL and VSR
549*
550 IF( ilvsl )
551 $ CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
552 $ rwork( iright ), n, vsl, ldvsl, ierr )
553 IF( ilvsr )
554 $ CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
555 $ rwork( iright ), n, vsr, ldvsr, ierr )
556*
557* Undo scaling
558*
559 IF( ilascl ) THEN
560 CALL clascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
561 CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
562 END IF
563*
564 IF( ilbscl ) THEN
565 CALL clascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
566 CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
567 END IF
568*
569 IF( wantst ) THEN
570*
571* Check if reordering is correct
572*
573 lastsl = .true.
574 sdim = 0
575 DO 20 i = 1, n
576 cursl = selctg( alpha( i ), beta( i ) )
577 IF( cursl )
578 $ sdim = sdim + 1
579 IF( cursl .AND. .NOT.lastsl )
580 $ info = n + 2
581 lastsl = cursl
582 20 CONTINUE
583*
584 END IF
585*
586 30 CONTINUE
587*
588 work( 1 ) = cmplx( lwkopt )
589*
590 RETURN
591*
592* End of CGGES3
593*
claqz0
recursive subroutine claqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, rec, info)
CLAQZ0
Definition claqz0.f:284
cgghd3
subroutine cgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
CGGHD3
Definition cgghd3.f:231

◆ cggesx()

subroutine cggesx ( character jobvsl,
character jobvsr,
character sort,
external selctg,
character sense,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
integer sdim,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvsl, * ) vsl,
integer ldvsl,
complex, dimension( ldvsr, * ) vsr,
integer ldvsr,
real, dimension( 2 ) rconde,
real, dimension( 2 ) rcondv,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer liwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> CGGESX 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 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 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 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 array, dimension (N)
!>
[out]BETA
!>          BETA is COMPLEX 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 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 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 REAL 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 REAL 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 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 REAL 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 WORK.
!>          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 CHGEQZ
!>                =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 CTGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 326 of file cggesx.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 REAL RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
344 COMPLEX 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 REAL ZERO, ONE
357 parameter( zero = 0.0e+0, one = 1.0e+0 )
358 COMPLEX CZERO, CONE
359 parameter( czero = ( 0.0e+0, 0.0e+0 ),
360 $ cone = ( 1.0e+0, 0.0e+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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
369 $ PR, SMLNUM
370* ..
371* .. Local Arrays ..
372 REAL DIF( 2 )
373* ..
374* .. External Subroutines ..
375 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
376 $ clascl, claset, ctgsen, cungqr, cunmqr, slabad,
377 $ xerbla
378* ..
379* .. External Functions ..
380 LOGICAL LSAME
381 INTEGER ILAENV
382 REAL CLANGE, SLAMCH
383 EXTERNAL lsame, ilaenv, clange, slamch
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, 'CGEQRF', ' ', n, 1, n, 0 ) )
465 maxwrk = max( maxwrk, n*( 1 +
466 $ ilaenv( 1, 'CUNMQR', ' ', n, 1, n, -1 ) ) )
467 IF( ilvsl ) THEN
468 maxwrk = max( maxwrk, n*( 1 +
469 $ ilaenv( 1, 'CUNGQR', ' ', 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( 'CGGESX', -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 = slamch( 'P' )
511 smlnum = slamch( 'S' )
512 bignum = one / smlnum
513 CALL slabad( 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 = clange( '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 clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
530*
531* Scale B if max element outside range [SMLNUM,BIGNUM]
532*
533 bnrm = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vsl, ldvsl )
576 IF( irows.GT.1 ) THEN
577 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
578 $ vsl( ilo+1, ilo ), ldvsl )
579 END IF
580 CALL cungqr( 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 claset( 'Full', n, n, czero, cone, vsr, ldvsr )
588*
589* Reduce to generalized Hessenberg form
590* (Workspace: none needed)
591*
592 CALL cgghrd( 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 chgeqz( '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 clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
625 IF( ilbscl )
626 $ CALL clascl( '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 ctgsen( 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 cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
671 $ rwork( iright ), n, vsl, ldvsl, ierr )
672*
673 IF( ilvsr )
674 $ CALL cggbak( '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 clascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
681 CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
682 END IF
683*
684 IF( ilbscl ) THEN
685 CALL clascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
686 CALL clascl( '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 CGGESX
714*

◆ cggev()

subroutine cggev ( character jobvl,
character jobvr,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvl, * ) vl,
integer ldvl,
complex, dimension( ldvr, * ) vr,
integer ldvr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CGGEV 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 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 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 array, dimension (N)
!>
[out]BETA
!>          BETA is COMPLEX 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 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 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 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 REAL 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 CHGEQZ,
!>                =N+2: error return from CTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 215 of file cggev.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 REAL RWORK( * )
228 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
229 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
230 $ WORK( * )
231* ..
232*
233* =====================================================================
234*
235* .. Parameters ..
236 REAL ZERO, ONE
237 parameter( zero = 0.0e0, one = 1.0e0 )
238 COMPLEX CZERO, CONE
239 parameter( czero = ( 0.0e0, 0.0e0 ),
240 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
249 $ SMLNUM, TEMP
250 COMPLEX X
251* ..
252* .. Local Arrays ..
253 LOGICAL LDUMMA( 1 )
254* ..
255* .. External Subroutines ..
256 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
257 $ clascl, claset, ctgevc, cungqr, cunmqr, slabad,
258 $ xerbla
259* ..
260* .. External Functions ..
261 LOGICAL LSAME
262 INTEGER ILAENV
263 REAL CLANGE, SLAMCH
264 EXTERNAL lsame, ilaenv, clange, slamch
265* ..
266* .. Intrinsic Functions ..
267 INTRINSIC abs, aimag, max, real, sqrt
268* ..
269* .. Statement Functions ..
270 REAL ABS1
271* ..
272* .. Statement Function definitions ..
273 abs1( x ) = abs( real( x ) ) + abs( aimag( 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, 'CGEQRF', ' ', n, 1, n, 0 ) )
333 lwkopt = max( lwkopt, n +
334 $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
335 IF( ilvl ) THEN
336 lwkopt = max( lwkopt, n +
337 $ n*ilaenv( 1, 'CUNGQR', ' ', 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( 'CGGEV ', -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 = slamch( 'E' )*slamch( 'B' )
360 smlnum = slamch( 'S' )
361 bignum = one / smlnum
362 CALL slabad( 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 = clange( '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 clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
379*
380* Scale B if max element outside range [SMLNUM,BIGNUM]
381*
382 bnrm = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vl, ldvl )
429 IF( irows.GT.1 ) THEN
430 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
431 $ vl( ilo+1, ilo ), ldvl )
432 END IF
433 CALL cungqr( 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 claset( '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 cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
449 $ ldvl, vr, ldvr, ierr )
450 ELSE
451 CALL cgghrd( '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 chgeqz( 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 ctgevc( 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 cggbak( '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 cggbak( '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 clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
546*
547 IF( ilbscl )
548 $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
549*
550 work( 1 ) = lwkopt
551 RETURN
552*
553* End of CGGEV
554*

◆ cggev3()

subroutine cggev3 ( character jobvl,
character jobvr,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvl, * ) vl,
integer ldvl,
complex, dimension( ldvr, * ) vr,
integer ldvr,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CGGEV3 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 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 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 array, dimension (N)
!>
[out]BETA
!>          BETA is COMPLEX 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 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 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 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 REAL 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 CHGEQZ,
!>                =N+2: error return from CTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 214 of file cggev3.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 REAL RWORK( * )
227 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
228 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
229 $ WORK( * )
230* ..
231*
232* =====================================================================
233*
234* .. Parameters ..
235 REAL ZERO, ONE
236 parameter( zero = 0.0e0, one = 1.0e0 )
237 COMPLEX CZERO, CONE
238 parameter( czero = ( 0.0e0, 0.0e0 ),
239 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
248 $ SMLNUM, TEMP
249 COMPLEX X
250* ..
251* .. Local Arrays ..
252 LOGICAL LDUMMA( 1 )
253* ..
254* .. External Subroutines ..
255 EXTERNAL cgeqrf, cggbak, cggbal, cgghd3, claqz0, clacpy,
256 $ clascl, claset, ctgevc, cungqr, cunmqr, slabad,
257 $ xerbla
258* ..
259* .. External Functions ..
260 LOGICAL LSAME
261 REAL CLANGE, SLAMCH
262 EXTERNAL lsame, clange, slamch
263* ..
264* .. Intrinsic Functions ..
265 INTRINSIC abs, aimag, max, real, sqrt
266* ..
267* .. Statement Functions ..
268 REAL ABS1
269* ..
270* .. Statement Function definitions ..
271 abs1( x ) = abs( real( x ) ) + abs( aimag( 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 cgeqrf( n, n, b, ldb, work, work, -1, ierr )
326 lwkopt = max( n, n+int( work( 1 ) ) )
327 CALL cunmqr( '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 cungqr( 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 cgghd3( 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 claqz0( '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 cgghd3( 'N', 'N', n, 1, n, a, lda, b, ldb, vl, ldvl,
344 $ vr, ldvr, work, -1, ierr )
345 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
346 CALL claqz0( '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 ) = cmplx( lwkopt )
352 END IF
353*
354 IF( info.NE.0 ) THEN
355 CALL xerbla( 'CGGEV3 ', -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 = slamch( 'E' )*slamch( 'B' )
369 smlnum = slamch( 'S' )
370 bignum = one / smlnum
371 CALL slabad( 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 = clange( '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 clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
388*
389* Scale B if max element outside range [SMLNUM,BIGNUM]
390*
391 bnrm = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vl, ldvl )
434 IF( irows.GT.1 ) THEN
435 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
436 $ vl( ilo+1, ilo ), ldvl )
437 END IF
438 CALL cungqr( 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 claset( '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 cgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
454 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk,
455 $ ierr )
456 ELSE
457 CALL cgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
458 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
459 $ work( iwrk ), lwork+1-iwrk, ierr )
460 END IF
461*
462* Perform QZ algorithm (Compute eigenvalues, and optionally, the
463* Schur form and Schur vectors)
464*
465 iwrk = itau
466 IF( ilv ) THEN
467 chtemp = 'S'
468 ELSE
469 chtemp = 'E'
470 END IF
471 CALL claqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
472 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
473 $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
474 IF( ierr.NE.0 ) THEN
475 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
476 info = ierr
477 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
478 info = ierr - n
479 ELSE
480 info = n + 1
481 END IF
482 GO TO 70
483 END IF
484*
485* Compute Eigenvectors
486*
487 IF( ilv ) THEN
488 IF( ilvl ) THEN
489 IF( ilvr ) THEN
490 chtemp = 'B'
491 ELSE
492 chtemp = 'L'
493 END IF
494 ELSE
495 chtemp = 'R'
496 END IF
497*
498 CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
499 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
500 $ ierr )
501 IF( ierr.NE.0 ) THEN
502 info = n + 2
503 GO TO 70
504 END IF
505*
506* Undo balancing on VL and VR and normalization
507*
508 IF( ilvl ) THEN
509 CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
510 $ rwork( iright ), n, vl, ldvl, ierr )
511 DO 30 jc = 1, n
512 temp = zero
513 DO 10 jr = 1, n
514 temp = max( temp, abs1( vl( jr, jc ) ) )
515 10 CONTINUE
516 IF( temp.LT.smlnum )
517 $ GO TO 30
518 temp = one / temp
519 DO 20 jr = 1, n
520 vl( jr, jc ) = vl( jr, jc )*temp
521 20 CONTINUE
522 30 CONTINUE
523 END IF
524 IF( ilvr ) THEN
525 CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
526 $ rwork( iright ), n, vr, ldvr, ierr )
527 DO 60 jc = 1, n
528 temp = zero
529 DO 40 jr = 1, n
530 temp = max( temp, abs1( vr( jr, jc ) ) )
531 40 CONTINUE
532 IF( temp.LT.smlnum )
533 $ GO TO 60
534 temp = one / temp
535 DO 50 jr = 1, n
536 vr( jr, jc ) = vr( jr, jc )*temp
537 50 CONTINUE
538 60 CONTINUE
539 END IF
540 END IF
541*
542* Undo scaling if necessary
543*
544 70 CONTINUE
545*
546 IF( ilascl )
547 $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
548*
549 IF( ilbscl )
550 $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
551*
552 work( 1 ) = cmplx( lwkopt )
553 RETURN
554*
555* End of CGGEV3
556*

◆ cggevx()

subroutine cggevx ( character balanc,
character jobvl,
character jobvr,
character sense,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( * ) alpha,
complex, dimension( * ) beta,
complex, dimension( ldvl, * ) vl,
integer ldvl,
complex, dimension( ldvr, * ) vr,
integer ldvr,
integer ilo,
integer ihi,
real, dimension( * ) lscale,
real, dimension( * ) rscale,
real abnrm,
real bbnrm,
real, dimension( * ) rconde,
real, dimension( * ) rcondv,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer, dimension( * ) iwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> CGGEVX 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 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 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 array, dimension (N)
!>
[out]BETA
!>          BETA is COMPLEX 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 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 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 REAL 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 REAL 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 REAL
!>          The one-norm of the balanced matrix A.
!>
[out]BBNRM
!>          BBNRM is REAL
!>          The one-norm of the balanced matrix B.
!>
[out]RCONDE
!>          RCONDE is REAL 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 REAL array, dimension (N)
!>          If SENSE = '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 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 REAL 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 CHGEQZ.
!>                =N+2: error return from CTGEVC.
!>
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 cggevx.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 REAL ABNRM, BBNRM
383* ..
384* .. Array Arguments ..
385 LOGICAL BWORK( * )
386 INTEGER IWORK( * )
387 REAL LSCALE( * ), RCONDE( * ), RCONDV( * ),
388 $ RSCALE( * ), RWORK( * )
389 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
390 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
391 $ WORK( * )
392* ..
393*
394* =====================================================================
395*
396* .. Parameters ..
397 REAL ZERO, ONE
398 parameter( zero = 0.0e+0, one = 1.0e+0 )
399 COMPLEX CZERO, CONE
400 parameter( czero = ( 0.0e+0, 0.0e+0 ),
401 $ cone = ( 1.0e+0, 0.0e+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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
410 $ SMLNUM, TEMP
411 COMPLEX X
412* ..
413* .. Local Arrays ..
414 LOGICAL LDUMMA( 1 )
415* ..
416* .. External Subroutines ..
417 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
418 $ clascl, claset, ctgevc, ctgsna, cungqr, cunmqr,
419 $ slabad, slascl, xerbla
420* ..
421* .. External Functions ..
422 LOGICAL LSAME
423 INTEGER ILAENV
424 REAL CLANGE, SLAMCH
425 EXTERNAL lsame, ilaenv, clange, slamch
426* ..
427* .. Intrinsic Functions ..
428 INTRINSIC abs, aimag, max, real, sqrt
429* ..
430* .. Statement Functions ..
431 REAL ABS1
432* ..
433* .. Statement Function definitions ..
434 abs1( x ) = abs( real( x ) ) + abs( aimag( 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, 'CGEQRF', ' ', n, 1, n, 0 ) )
517 maxwrk = max( maxwrk,
518 $ n + n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
519 IF( ilvl ) THEN
520 maxwrk = max( maxwrk, n +
521 $ n*ilaenv( 1, 'CUNGQR', ' ', 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( 'CGGEVX', -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 = slamch( 'P' )
546 smlnum = slamch( 'S' )
547 bignum = one / smlnum
548 CALL slabad( 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 = clange( '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 clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
565*
566* Scale B if max element outside range [SMLNUM,BIGNUM]
567*
568 bnrm = clange( '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 clascl( '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 cggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
584 $ rwork, ierr )
585*
586* Compute ABNRM and BBNRM
587*
588 abnrm = clange( '1', n, n, a, lda, rwork( 1 ) )
589 IF( ilascl ) THEN
590 rwork( 1 ) = abnrm
591 CALL slascl( 'G', 0, 0, anrmto, anrm, 1, 1, rwork( 1 ), 1,
592 $ ierr )
593 abnrm = rwork( 1 )
594 END IF
595*
596 bbnrm = clange( '1', n, n, b, ldb, rwork( 1 ) )
597 IF( ilbscl ) THEN
598 rwork( 1 ) = bbnrm
599 CALL slascl( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vl, ldvl )
630 IF( irows.GT.1 ) THEN
631 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
632 $ vl( ilo+1, ilo ), ldvl )
633 END IF
634 CALL cungqr( 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 claset( '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 cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
649 $ ldvl, vr, ldvr, ierr )
650 ELSE
651 CALL cgghrd( '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 chgeqz( 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* CTGEVC: (Complex Workspace: need 2*N )
683* (Real Workspace: need 2*N )
684* CTGSNA: (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 ctgevc( 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 (CTGEVC) and estimate condition
711* numbers (CTGSNA). 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 ctgevc( '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 ctgsna( 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 cggbak( 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 cggbak( 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 clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
792*
793 IF( ilbscl )
794 $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
795*
796 work( 1 ) = maxwrk
797 RETURN
798*
799* End of CGGEVX
800*
ctgsna
subroutine ctgsna(job, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, ldvr, s, dif, mm, m, work, lwork, iwork, info)
CTGSNA
Definition ctgsna.f:311