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

Functions

subroutine dgegs (jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, info)
  DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine dgegv (jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
  DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine dgees (jobvs, sort, select, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info)
  DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine dgeesx (jobvs, sort, select, sense, n, a, lda, sdim, wr, wi, vs, ldvs, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
  DGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine dgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
  DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine dgeevx (balanc, jobvl, jobvr, sense, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv, work, lwork, iwork, info)
  DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine dgges (jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
  DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine dgges3 (jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
  DGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)
subroutine dggesx (jobvsl, jobvsr, sort, selctg, sense, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
  DGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices
subroutine dggev (jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
  DGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine dggev3 (jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
  DGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)
subroutine dggevx (balanc, jobvl, jobvr, sense, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)
  DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

Detailed Description

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

Function Documentation

◆ dgees()

subroutine dgees ( character jobvs,
character sort,
external select,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
integer sdim,
double precision, dimension( * ) wr,
double precision, dimension( * ) wi,
double precision, dimension( ldvs, * ) vs,
integer ldvs,
double precision, dimension( * ) work,
integer lwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> DGEES computes for an N-by-N real nonsymmetric matrix A, the
!> eigenvalues, the real Schur form T, and, optionally, the matrix of
!> Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
!>
!> Optionally, it also orders the eigenvalues on the diagonal of the
!> real Schur form so that selected eigenvalues are at the top left.
!> The leading columns of Z then form an orthonormal basis for the
!> invariant subspace corresponding to the selected eigenvalues.
!>
!> A matrix is in real Schur form if it is upper quasi-triangular with
!> 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
!> form
!>         [  a  b  ]
!>         [  c  a  ]
!>
!> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
!>
Parameters
[in]JOBVS
!>          JOBVS is CHARACTER*1
!>          = 'N': Schur vectors are not computed;
!>          = 'V': Schur vectors are computed.
!>
[in]SORT
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the Schur form.
!>          = 'N': Eigenvalues are not ordered;
!>          = 'S': Eigenvalues are ordered (see SELECT).
!>
[in]SELECT
!>          SELECT is a LOGICAL FUNCTION of two DOUBLE PRECISION arguments
!>          SELECT must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'S', SELECT is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          If SORT = 'N', SELECT is not referenced.
!>          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
!>          SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
!>          conjugate pair of eigenvalues is selected, then both complex
!>          eigenvalues are selected.
!>          Note that a selected complex eigenvalue may no longer
!>          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
!>          ordering may change the value of complex eigenvalues
!>          (especially if the eigenvalue is ill-conditioned); in this
!>          case INFO is set to N+2 (see INFO below).
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A. N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the N-by-N matrix A.
!>          On exit, A has been overwritten by its real Schur form T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]SDIM
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>                         for which SELECT is true. (Complex conjugate
!>                         pairs for which SELECT is true for either
!>                         eigenvalue count as 2.)
!>
[out]WR
!>          WR is DOUBLE PRECISION array, dimension (N)
!>
[out]WI
!>          WI is DOUBLE PRECISION array, dimension (N)
!>          WR and WI contain the real and imaginary parts,
!>          respectively, of the computed eigenvalues in the same order
!>          that they appear on the diagonal of the output Schur form T.
!>          Complex conjugate pairs of eigenvalues will appear
!>          consecutively with the eigenvalue having the positive
!>          imaginary part first.
!>
[out]VS
!>          VS is DOUBLE PRECISION array, dimension (LDVS,N)
!>          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
!>          vectors.
!>          If JOBVS = 'N', VS is not referenced.
!>
[in]LDVS
!>          LDVS is INTEGER
!>          The leading dimension of the array VS.  LDVS >= 1; if
!>          JOBVS = 'V', LDVS >= N.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,3*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]BWORK
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value.
!>          > 0: if INFO = i, and i is
!>             <= N: the QR algorithm failed to compute all the
!>                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
!>                   contain those eigenvalues which have converged; if
!>                   JOBVS = 'V', VS contains the 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 214 of file dgees.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 JOBVS, SORT
223 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
224* ..
225* .. Array Arguments ..
226 LOGICAL BWORK( * )
227 DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
228 $ WR( * )
229* ..
230* .. Function Arguments ..
231 LOGICAL SELECT
232 EXTERNAL SELECT
233* ..
234*
235* =====================================================================
236*
237* .. Parameters ..
238 DOUBLE PRECISION ZERO, ONE
239 parameter( zero = 0.0d0, one = 1.0d0 )
240* ..
241* .. Local Scalars ..
242 LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
243 $ WANTVS
244 INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
245 $ IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
246 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
247* ..
248* .. Local Arrays ..
249 INTEGER IDUM( 1 )
250 DOUBLE PRECISION DUM( 1 )
251* ..
252* .. External Subroutines ..
253 EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
254 $ dlabad, dlascl, dorghr, dswap, dtrsen, xerbla
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 INTEGER ILAENV
259 DOUBLE PRECISION DLAMCH, DLANGE
260 EXTERNAL lsame, ilaenv, dlamch, dlange
261* ..
262* .. Intrinsic Functions ..
263 INTRINSIC max, sqrt
264* ..
265* .. Executable Statements ..
266*
267* Test the input arguments
268*
269 info = 0
270 lquery = ( lwork.EQ.-1 )
271 wantvs = lsame( jobvs, 'V' )
272 wantst = lsame( sort, 'S' )
273 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
274 info = -1
275 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
276 info = -2
277 ELSE IF( n.LT.0 ) THEN
278 info = -4
279 ELSE IF( lda.LT.max( 1, n ) ) THEN
280 info = -6
281 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
282 info = -11
283 END IF
284*
285* Compute workspace
286* (Note: Comments in the code beginning "Workspace:" describe the
287* minimal amount of workspace needed at that point in the code,
288* as well as the preferred amount for good performance.
289* NB refers to the optimal block size for the immediately
290* following subroutine, as returned by ILAENV.
291* HSWORK refers to the workspace preferred by DHSEQR, as
292* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
293* the worst case.)
294*
295 IF( info.EQ.0 ) THEN
296 IF( n.EQ.0 ) THEN
297 minwrk = 1
298 maxwrk = 1
299 ELSE
300 maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
301 minwrk = 3*n
302*
303 CALL dhseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
304 $ work, -1, ieval )
305 hswork = work( 1 )
306*
307 IF( .NOT.wantvs ) THEN
308 maxwrk = max( maxwrk, n + hswork )
309 ELSE
310 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
311 $ 'DORGHR', ' ', n, 1, n, -1 ) )
312 maxwrk = max( maxwrk, n + hswork )
313 END IF
314 END IF
315 work( 1 ) = maxwrk
316*
317 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
318 info = -13
319 END IF
320 END IF
321*
322 IF( info.NE.0 ) THEN
323 CALL xerbla( 'DGEES ', -info )
324 RETURN
325 ELSE IF( lquery ) THEN
326 RETURN
327 END IF
328*
329* Quick return if possible
330*
331 IF( n.EQ.0 ) THEN
332 sdim = 0
333 RETURN
334 END IF
335*
336* Get machine constants
337*
338 eps = dlamch( 'P' )
339 smlnum = dlamch( 'S' )
340 bignum = one / smlnum
341 CALL dlabad( smlnum, bignum )
342 smlnum = sqrt( smlnum ) / eps
343 bignum = one / smlnum
344*
345* Scale A if max element outside range [SMLNUM,BIGNUM]
346*
347 anrm = dlange( 'M', n, n, a, lda, dum )
348 scalea = .false.
349 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
350 scalea = .true.
351 cscale = smlnum
352 ELSE IF( anrm.GT.bignum ) THEN
353 scalea = .true.
354 cscale = bignum
355 END IF
356 IF( scalea )
357 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
358*
359* Permute the matrix to make it more nearly triangular
360* (Workspace: need N)
361*
362 ibal = 1
363 CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
364*
365* Reduce to upper Hessenberg form
366* (Workspace: need 3*N, prefer 2*N+N*NB)
367*
368 itau = n + ibal
369 iwrk = n + itau
370 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
371 $ lwork-iwrk+1, ierr )
372*
373 IF( wantvs ) THEN
374*
375* Copy Householder vectors to VS
376*
377 CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
378*
379* Generate orthogonal matrix in VS
380* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
381*
382 CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
383 $ lwork-iwrk+1, ierr )
384 END IF
385*
386 sdim = 0
387*
388* Perform QR iteration, accumulating Schur vectors in VS if desired
389* (Workspace: need N+1, prefer N+HSWORK (see comments) )
390*
391 iwrk = itau
392 CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
393 $ work( iwrk ), lwork-iwrk+1, ieval )
394 IF( ieval.GT.0 )
395 $ info = ieval
396*
397* Sort eigenvalues if desired
398*
399 IF( wantst .AND. info.EQ.0 ) THEN
400 IF( scalea ) THEN
401 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
402 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
403 END IF
404 DO 10 i = 1, n
405 bwork( i ) = SELECT( wr( i ), wi( i ) )
406 10 CONTINUE
407*
408* Reorder eigenvalues and transform Schur vectors
409* (Workspace: none needed)
410*
411 CALL dtrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
412 $ sdim, s, sep, work( iwrk ), lwork-iwrk+1, idum, 1,
413 $ icond )
414 IF( icond.GT.0 )
415 $ info = n + icond
416 END IF
417*
418 IF( wantvs ) THEN
419*
420* Undo balancing
421* (Workspace: need N)
422*
423 CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
424 $ ierr )
425 END IF
426*
427 IF( scalea ) THEN
428*
429* Undo scaling for the Schur form of A
430*
431 CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
432 CALL dcopy( n, a, lda+1, wr, 1 )
433 IF( cscale.EQ.smlnum ) THEN
434*
435* If scaling back towards underflow, adjust WI if an
436* offdiagonal element of a 2-by-2 block in the Schur form
437* underflows.
438*
439 IF( ieval.GT.0 ) THEN
440 i1 = ieval + 1
441 i2 = ihi - 1
442 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi,
443 $ max( ilo-1, 1 ), ierr )
444 ELSE IF( wantst ) THEN
445 i1 = 1
446 i2 = n - 1
447 ELSE
448 i1 = ilo
449 i2 = ihi - 1
450 END IF
451 inxt = i1 - 1
452 DO 20 i = i1, i2
453 IF( i.LT.inxt )
454 $ GO TO 20
455 IF( wi( i ).EQ.zero ) THEN
456 inxt = i + 1
457 ELSE
458 IF( a( i+1, i ).EQ.zero ) THEN
459 wi( i ) = zero
460 wi( i+1 ) = zero
461 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
462 $ zero ) THEN
463 wi( i ) = zero
464 wi( i+1 ) = zero
465 IF( i.GT.1 )
466 $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
467 IF( n.GT.i+1 )
468 $ CALL dswap( n-i-1, a( i, i+2 ), lda,
469 $ a( i+1, i+2 ), lda )
470 IF( wantvs ) THEN
471 CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
472 END IF
473 a( i, i+1 ) = a( i+1, i )
474 a( i+1, i ) = zero
475 END IF
476 inxt = i + 2
477 END IF
478 20 CONTINUE
479 END IF
480*
481* Undo scaling for the imaginary part of the eigenvalues
482*
483 CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
484 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
485 END IF
486*
487 IF( wantst .AND. info.EQ.0 ) THEN
488*
489* Check if reordering successful
490*
491 lastsl = .true.
492 lst2sl = .true.
493 sdim = 0
494 ip = 0
495 DO 30 i = 1, n
496 cursl = SELECT( wr( i ), wi( i ) )
497 IF( wi( i ).EQ.zero ) THEN
498 IF( cursl )
499 $ sdim = sdim + 1
500 ip = 0
501 IF( cursl .AND. .NOT.lastsl )
502 $ info = n + 2
503 ELSE
504 IF( ip.EQ.1 ) THEN
505*
506* Last eigenvalue of conjugate pair
507*
508 cursl = cursl .OR. lastsl
509 lastsl = cursl
510 IF( cursl )
511 $ sdim = sdim + 2
512 ip = -1
513 IF( cursl .AND. .NOT.lst2sl )
514 $ info = n + 2
515 ELSE
516*
517* First eigenvalue of conjugate pair
518*
519 ip = 1
520 END IF
521 END IF
522 lst2sl = lastsl
523 lastsl = cursl
524 30 CONTINUE
525 END IF
526*
527 work( 1 ) = maxwrk
528 RETURN
529*
530* End of DGEES
531*
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:114
dgehrd
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
dgebal
subroutine dgebal(job, n, a, lda, ilo, ihi, scale, info)
DGEBAL
Definition dgebal.f:160
dgebak
subroutine dgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
DGEBAK
Definition dgebak.f:130
dhseqr
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:316
dtrsen
subroutine dtrsen(job, compq, select, n, t, ldt, q, ldq, wr, wi, m, s, sep, work, lwork, iwork, liwork, info)
DTRSEN
Definition dtrsen.f:313
dorghr
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:126
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21

◆ dgeesx()

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

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

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

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

Definition at line 278 of file dgeesx.f.

281*
282* -- LAPACK driver routine --
283* -- LAPACK is a software package provided by Univ. of Tennessee, --
284* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
285*
286* .. Scalar Arguments ..
287 CHARACTER JOBVS, SENSE, SORT
288 INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
289 DOUBLE PRECISION RCONDE, RCONDV
290* ..
291* .. Array Arguments ..
292 LOGICAL BWORK( * )
293 INTEGER IWORK( * )
294 DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
295 $ WR( * )
296* ..
297* .. Function Arguments ..
298 LOGICAL SELECT
299 EXTERNAL SELECT
300* ..
301*
302* =====================================================================
303*
304* .. Parameters ..
305 DOUBLE PRECISION ZERO, ONE
306 parameter( zero = 0.0d0, one = 1.0d0 )
307* ..
308* .. Local Scalars ..
309 LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
310 $ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
311 INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
312 $ IHI, ILO, INXT, IP, ITAU, IWRK, LIWRK, LWRK,
313 $ MAXWRK, MINWRK
314 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
315* ..
316* .. Local Arrays ..
317 DOUBLE PRECISION DUM( 1 )
318* ..
319* .. External Subroutines ..
320 EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
321 $ dlascl, dorghr, dswap, dtrsen, xerbla
322* ..
323* .. External Functions ..
324 LOGICAL LSAME
325 INTEGER ILAENV
326 DOUBLE PRECISION DLAMCH, DLANGE
327 EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
328* ..
329* .. Intrinsic Functions ..
330 INTRINSIC max, sqrt
331* ..
332* .. Executable Statements ..
333*
334* Test the input arguments
335*
336 info = 0
337 wantvs = lsame( jobvs, 'V' )
338 wantst = lsame( sort, 'S' )
339 wantsn = lsame( sense, 'N' )
340 wantse = lsame( sense, 'E' )
341 wantsv = lsame( sense, 'V' )
342 wantsb = lsame( sense, 'B' )
343 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
344*
345 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
346 info = -1
347 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
348 info = -2
349 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
350 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
351 info = -4
352 ELSE IF( n.LT.0 ) THEN
353 info = -5
354 ELSE IF( lda.LT.max( 1, n ) ) THEN
355 info = -7
356 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
357 info = -12
358 END IF
359*
360* Compute workspace
361* (Note: Comments in the code beginning "RWorkspace:" describe the
362* minimal amount of real workspace needed at that point in the
363* code, as well as the preferred amount for good performance.
364* IWorkspace refers to integer workspace.
365* NB refers to the optimal block size for the immediately
366* following subroutine, as returned by ILAENV.
367* HSWORK refers to the workspace preferred by DHSEQR, as
368* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
369* the worst case.
370* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
371* depends on SDIM, which is computed by the routine DTRSEN later
372* in the code.)
373*
374 IF( info.EQ.0 ) THEN
375 liwrk = 1
376 IF( n.EQ.0 ) THEN
377 minwrk = 1
378 lwrk = 1
379 ELSE
380 maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
381 minwrk = 3*n
382*
383 CALL dhseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
384 $ work, -1, ieval )
385 hswork = work( 1 )
386*
387 IF( .NOT.wantvs ) THEN
388 maxwrk = max( maxwrk, n + hswork )
389 ELSE
390 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
391 $ 'DORGHR', ' ', n, 1, n, -1 ) )
392 maxwrk = max( maxwrk, n + hswork )
393 END IF
394 lwrk = maxwrk
395 IF( .NOT.wantsn )
396 $ lwrk = max( lwrk, n + ( n*n )/2 )
397 IF( wantsv .OR. wantsb )
398 $ liwrk = ( n*n )/4
399 END IF
400 iwork( 1 ) = liwrk
401 work( 1 ) = lwrk
402*
403 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
404 info = -16
405 ELSE IF( liwork.LT.1 .AND. .NOT.lquery ) THEN
406 info = -18
407 END IF
408 END IF
409*
410 IF( info.NE.0 ) THEN
411 CALL xerbla( 'DGEESX', -info )
412 RETURN
413 ELSE IF( lquery ) THEN
414 RETURN
415 END IF
416*
417* Quick return if possible
418*
419 IF( n.EQ.0 ) THEN
420 sdim = 0
421 RETURN
422 END IF
423*
424* Get machine constants
425*
426 eps = dlamch( 'P' )
427 smlnum = dlamch( 'S' )
428 bignum = one / smlnum
429 CALL dlabad( smlnum, bignum )
430 smlnum = sqrt( smlnum ) / eps
431 bignum = one / smlnum
432*
433* Scale A if max element outside range [SMLNUM,BIGNUM]
434*
435 anrm = dlange( 'M', n, n, a, lda, dum )
436 scalea = .false.
437 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
438 scalea = .true.
439 cscale = smlnum
440 ELSE IF( anrm.GT.bignum ) THEN
441 scalea = .true.
442 cscale = bignum
443 END IF
444 IF( scalea )
445 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
446*
447* Permute the matrix to make it more nearly triangular
448* (RWorkspace: need N)
449*
450 ibal = 1
451 CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
452*
453* Reduce to upper Hessenberg form
454* (RWorkspace: need 3*N, prefer 2*N+N*NB)
455*
456 itau = n + ibal
457 iwrk = n + itau
458 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
459 $ lwork-iwrk+1, ierr )
460*
461 IF( wantvs ) THEN
462*
463* Copy Householder vectors to VS
464*
465 CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
466*
467* Generate orthogonal matrix in VS
468* (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
469*
470 CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
471 $ lwork-iwrk+1, ierr )
472 END IF
473*
474 sdim = 0
475*
476* Perform QR iteration, accumulating Schur vectors in VS if desired
477* (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
478*
479 iwrk = itau
480 CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
481 $ work( iwrk ), lwork-iwrk+1, ieval )
482 IF( ieval.GT.0 )
483 $ info = ieval
484*
485* Sort eigenvalues if desired
486*
487 IF( wantst .AND. info.EQ.0 ) THEN
488 IF( scalea ) THEN
489 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
490 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
491 END IF
492 DO 10 i = 1, n
493 bwork( i ) = SELECT( wr( i ), wi( i ) )
494 10 CONTINUE
495*
496* Reorder eigenvalues, transform Schur vectors, and compute
497* reciprocal condition numbers
498* (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
499* otherwise, need N )
500* (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
501* otherwise, need 0 )
502*
503 CALL dtrsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
504 $ sdim, rconde, rcondv, work( iwrk ), lwork-iwrk+1,
505 $ iwork, liwork, icond )
506 IF( .NOT.wantsn )
507 $ maxwrk = max( maxwrk, n+2*sdim*( n-sdim ) )
508 IF( icond.EQ.-15 ) THEN
509*
510* Not enough real workspace
511*
512 info = -16
513 ELSE IF( icond.EQ.-17 ) THEN
514*
515* Not enough integer workspace
516*
517 info = -18
518 ELSE IF( icond.GT.0 ) THEN
519*
520* DTRSEN failed to reorder or to restore standard Schur form
521*
522 info = icond + n
523 END IF
524 END IF
525*
526 IF( wantvs ) THEN
527*
528* Undo balancing
529* (RWorkspace: need N)
530*
531 CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
532 $ ierr )
533 END IF
534*
535 IF( scalea ) THEN
536*
537* Undo scaling for the Schur form of A
538*
539 CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
540 CALL dcopy( n, a, lda+1, wr, 1 )
541 IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
542 dum( 1 ) = rcondv
543 CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
544 rcondv = dum( 1 )
545 END IF
546 IF( cscale.EQ.smlnum ) THEN
547*
548* If scaling back towards underflow, adjust WI if an
549* offdiagonal element of a 2-by-2 block in the Schur form
550* underflows.
551*
552 IF( ieval.GT.0 ) THEN
553 i1 = ieval + 1
554 i2 = ihi - 1
555 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
556 $ ierr )
557 ELSE IF( wantst ) THEN
558 i1 = 1
559 i2 = n - 1
560 ELSE
561 i1 = ilo
562 i2 = ihi - 1
563 END IF
564 inxt = i1 - 1
565 DO 20 i = i1, i2
566 IF( i.LT.inxt )
567 $ GO TO 20
568 IF( wi( i ).EQ.zero ) THEN
569 inxt = i + 1
570 ELSE
571 IF( a( i+1, i ).EQ.zero ) THEN
572 wi( i ) = zero
573 wi( i+1 ) = zero
574 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
575 $ zero ) THEN
576 wi( i ) = zero
577 wi( i+1 ) = zero
578 IF( i.GT.1 )
579 $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
580 IF( n.GT.i+1 )
581 $ CALL dswap( n-i-1, a( i, i+2 ), lda,
582 $ a( i+1, i+2 ), lda )
583 IF( wantvs ) THEN
584 CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
585 END IF
586 a( i, i+1 ) = a( i+1, i )
587 a( i+1, i ) = zero
588 END IF
589 inxt = i + 2
590 END IF
591 20 CONTINUE
592 END IF
593 CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
594 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
595 END IF
596*
597 IF( wantst .AND. info.EQ.0 ) THEN
598*
599* Check if reordering successful
600*
601 lastsl = .true.
602 lst2sl = .true.
603 sdim = 0
604 ip = 0
605 DO 30 i = 1, n
606 cursl = SELECT( wr( i ), wi( i ) )
607 IF( wi( i ).EQ.zero ) THEN
608 IF( cursl )
609 $ sdim = sdim + 1
610 ip = 0
611 IF( cursl .AND. .NOT.lastsl )
612 $ info = n + 2
613 ELSE
614 IF( ip.EQ.1 ) THEN
615*
616* Last eigenvalue of conjugate pair
617*
618 cursl = cursl .OR. lastsl
619 lastsl = cursl
620 IF( cursl )
621 $ sdim = sdim + 2
622 ip = -1
623 IF( cursl .AND. .NOT.lst2sl )
624 $ info = n + 2
625 ELSE
626*
627* First eigenvalue of conjugate pair
628*
629 ip = 1
630 END IF
631 END IF
632 lst2sl = lastsl
633 lastsl = cursl
634 30 CONTINUE
635 END IF
636*
637 work( 1 ) = maxwrk
638 IF( wantsv .OR. wantsb ) THEN
639 iwork( 1 ) = max( 1, sdim*( n-sdim ) )
640 ELSE
641 iwork( 1 ) = 1
642 END IF
643*
644 RETURN
645*
646* End of DGEESX
647*

◆ dgeev()

subroutine dgeev ( character jobvl,
character jobvr,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) wr,
double precision, dimension( * ) wi,
double precision, dimension( ldvl, * ) vl,
integer ldvl,
double precision, dimension( ldvr, * ) vr,
integer ldvr,
double precision, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> DGEEV computes for an N-by-N real 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 A 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 DOUBLE PRECISION 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]WR
!>          WR is DOUBLE PRECISION array, dimension (N)
!>
[out]WI
!>          WI is DOUBLE PRECISION array, dimension (N)
!>          WR and WI contain the real and imaginary parts,
!>          respectively, of the computed eigenvalues.  Complex
!>          conjugate pairs of eigenvalues appear consecutively
!>          with the eigenvalue having the positive imaginary part
!>          first.
!>
[out]VL
!>          VL is DOUBLE PRECISION 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.
!>          If the j-th eigenvalue is real, then u(j) = VL(:,j),
!>          the j-th column of VL.
!>          If the j-th and (j+1)-st eigenvalues form a complex
!>          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
!>          u(j+1) = VL(:,j) - i*VL(:,j+1).
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the array VL.  LDVL >= 1; if
!>          JOBVL = 'V', LDVL >= N.
!>
[out]VR
!>          VR is DOUBLE PRECISION 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.
!>          If the j-th eigenvalue is real, then v(j) = VR(:,j),
!>          the j-th column of VR.
!>          If the j-th and (j+1)-st eigenvalues form a complex
!>          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
!>          v(j+1) = VR(:,j) - i*VR(:,j+1).
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the array VR.  LDVR >= 1; if
!>          JOBVR = 'V', LDVR >= N.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,3*N), and
!>          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*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]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 WR and WI contain eigenvalues which
!>                have converged.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 190 of file dgeev.f.

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

◆ dgeevx()

subroutine dgeevx ( character balanc,
character jobvl,
character jobvr,
character sense,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) wr,
double precision, dimension( * ) wi,
double precision, dimension( ldvl, * ) vl,
integer ldvl,
double precision, dimension( ldvr, * ) vr,
integer ldvr,
integer ilo,
integer ihi,
double precision, dimension( * ) scale,
double precision abnrm,
double precision, dimension( * ) rconde,
double precision, dimension( * ) rcondv,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer info )

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

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

Purpose:
!>
!> DGEEVX computes for an N-by-N real 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, i.e. 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 DOUBLE PRECISION 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 real Schur form of the balanced
!>          version of the input matrix A.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]WR
!>          WR is DOUBLE PRECISION array, dimension (N)
!>
[out]WI
!>          WI is DOUBLE PRECISION array, dimension (N)
!>          WR and WI contain the real and imaginary parts,
!>          respectively, of the computed eigenvalues.  Complex
!>          conjugate pairs of eigenvalues will appear consecutively
!>          with the eigenvalue having the positive imaginary part
!>          first.
!>
[out]VL
!>          VL is DOUBLE PRECISION 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.
!>          If the j-th eigenvalue is real, then u(j) = VL(:,j),
!>          the j-th column of VL.
!>          If the j-th and (j+1)-st eigenvalues form a complex
!>          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
!>          u(j+1) = VL(:,j) - i*VL(:,j+1).
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of the array VL.  LDVL >= 1; if
!>          JOBVL = 'V', LDVL >= N.
!>
[out]VR
!>          VR is DOUBLE PRECISION 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.
!>          If the j-th eigenvalue is real, then v(j) = VR(:,j),
!>          the j-th column of VR.
!>          If the j-th and (j+1)-st eigenvalues form a complex
!>          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
!>          v(j+1) = VR(:,j) - i*VR(:,j+1).
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the array 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 determined when A was
!>          balanced.  The balanced A(i,j) = 0 if I > J and
!>          J = 1,...,ILO-1 or I = IHI+1,...,N.
!>
[out]SCALE
!>          SCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          when balancing A.  If P(j) is the index of the row and column
!>          interchanged with row and column j, and D(j) is the scaling
!>          factor applied to row and column j, then
!>          SCALE(J) = P(J),    for J = 1,...,ILO-1
!>                   = D(J),    for J = ILO,...,IHI
!>                   = P(J)     for J = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]ABNRM
!>          ABNRM is DOUBLE PRECISION
!>          The one-norm of the balanced matrix (the maximum
!>          of the sum of absolute values of elements of any column).
!>
[out]RCONDE
!>          RCONDE is DOUBLE PRECISION array, dimension (N)
!>          RCONDE(j) is the reciprocal condition number of the j-th
!>          eigenvalue.
!>
[out]RCONDV
!>          RCONDV is DOUBLE PRECISION array, dimension (N)
!>          RCONDV(j) is the reciprocal condition number of the j-th
!>          right eigenvector.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.   If SENSE = 'N' or 'E',
!>          LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
!>          LWORK >= 3*N.  If SENSE = 'V' or 'B', LWORK >= N*(N+6).
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (2*N-2)
!>          If SENSE = 'N' or 'E', not referenced.
!>
[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 WR
!>                and WI contain eigenvalues which have converged.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 303 of file dgeevx.f.

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

◆ dgegs()

subroutine dgegs ( character jobvsl,
character jobvsr,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvsl, * ) vsl,
integer ldvsl,
double precision, dimension( ldvsr, * ) vsr,
integer ldvsr,
double precision, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine DGGES.
!>
!> DGEGS computes the eigenvalues, real Schur form, and, optionally,
!> left and or/right Schur vectors of a real matrix pair (A,B).
!> Given two square matrices A and B, the generalized real Schur
!> factorization has the form
!>
!>   A = Q*S*Z**T,  B = Q*T*Z**T
!>
!> where Q and Z are orthogonal matrices, T is upper triangular, and S
!> is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
!> blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
!> of eigenvalues of (A,B).  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
!> DGEGV should be used instead.  See DGEGV 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 DOUBLE PRECISION array, dimension (LDA, N)
!>          On entry, the matrix A.
!>          On exit, the upper quasi-triangular matrix S from the
!>          generalized real Schur factorization.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB, N)
!>          On entry, the matrix B.
!>          On exit, the upper triangular matrix T from the generalized
!>          real Schur factorization.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>          The real parts of each scalar alpha defining an eigenvalue
!>          of GNEP.
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>          The imaginary parts of each scalar alpha defining an
!>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
!>          eigenvalue is real; if positive, then the j-th and (j+1)-st
!>          eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) = -ALPHAI(j).
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          The scalars beta that define the eigenvalues of GNEP.
!>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,4*N).
!>          For good performance, LWORK must generally be larger.
!>          To compute the optimal value of LWORK, call ILAENV to get
!>          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
!>          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
!>          The optimal LWORK is  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]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 ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  errors that usually indicate LAPACK problems:
!>                =N+1: error return from DGGBAL
!>                =N+2: error return from DGEQRF
!>                =N+3: error return from DORMQR
!>                =N+4: error return from DORGQR
!>                =N+5: error return from DGGHRD
!>                =N+6: error return from DHGEQZ (other than failed
!>                                                iteration)
!>                =N+7: error return from DGGBAK (computing VSL)
!>                =N+8: error return from DGGBAK (computing VSR)
!>                =N+9: error return from DLASCL (various places)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file dgegs.f.

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

◆ dgegv()

subroutine dgegv ( character jobvl,
character jobvr,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvl, * ) vl,
integer ldvl,
double precision, dimension( ldvr, * ) vr,
integer ldvr,
double precision, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> This routine is deprecated and has been replaced by routine DGGEV.
!>
!> DGEGV computes the eigenvalues and, optionally, the left and/or right
!> eigenvectors of a real 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 DOUBLE PRECISION array, dimension (LDA, N)
!>          On entry, the matrix A.
!>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
!>          contains the real 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
!>          blocks from the Schur form will be correct.  See DGGHRD and
!>          DHGEQZ for details.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION 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 those elements of
!>          B corresponding to the diagonal blocks from the Schur form of
!>          A will be correct.  See DGGHRD and DHGEQZ for details.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>          The real parts of each scalar alpha defining an eigenvalue of
!>          GNEP.
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>          The imaginary parts of each scalar alpha defining an
!>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
!>          eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) = -ALPHAI(j).
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          The scalars beta that define the eigenvalues of GNEP.
!>
!>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(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 DOUBLE PRECISION 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.
!>          If the j-th eigenvalue is real, then u(j) = VL(:,j).
!>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
!>          pair, then
!>             u(j) = VL(:,j) + i*VL(:,j+1)
!>          and
!>            u(j+1) = VL(:,j) - i*VL(:,j+1).
!>
!>          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 DOUBLE PRECISION 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.
!>          If the j-th eigenvalue is real, then x(j) = VR(:,j).
!>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
!>          pair, then
!>            x(j) = VR(:,j) + i*VR(:,j+1)
!>          and
!>            x(j+1) = VR(:,j) - i*VR(:,j+1).
!>
!>          Each eigenvector is scaled so that its largest component has
!>          abs(real part) + abs(imag. part) = 1, except for eigenvalues
!>          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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,8*N).
!>          For good performance, LWORK must generally be larger.
!>          To compute the optimal value of LWORK, call ILAENV to get
!>          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
!>          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
!>          The optimal LWORK is:
!>              2*N + MAX( 6*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]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 ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  errors that usually indicate LAPACK problems:
!>                =N+1: error return from DGGBAL
!>                =N+2: error return from DGEQRF
!>                =N+3: error return from DORMQR
!>                =N+4: error return from DORGQR
!>                =N+5: error return from DGGHRD
!>                =N+6: error return from DHGEQZ (other than failed
!>                                                iteration)
!>                =N+7: error return from DTGEVC
!>                =N+8: error return from DGGBAK (computing VL)
!>                =N+9: error return from DGGBAK (computing VR)
!>                =N+10: error return from DLASCL (various calls)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Balancing
!>  ---------
!>
!>  This driver calls DGGBAL 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, DGGBAK 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 real Schur
!>  form[*] of the  versions of A and B.  If no eigenvectors
!>  are computed, then only the diagonal blocks will be correct.
!>
!>  [*] See DHGEQZ, DGEGS, or read the book ,
!>      by Golub & van Loan, pub. by Johns Hopkins U. Press.
!>

Definition at line 304 of file dgegv.f.

306*
307* -- LAPACK driver routine --
308* -- LAPACK is a software package provided by Univ. of Tennessee, --
309* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310*
311* .. Scalar Arguments ..
312 CHARACTER JOBVL, JOBVR
313 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
314* ..
315* .. Array Arguments ..
316 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
317 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
318 $ VR( LDVR, * ), WORK( * )
319* ..
320*
321* =====================================================================
322*
323* .. Parameters ..
324 DOUBLE PRECISION ZERO, ONE
325 parameter( zero = 0.0d0, one = 1.0d0 )
326* ..
327* .. Local Scalars ..
328 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329 CHARACTER CHTEMP
330 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
331 $ IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
332 $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
333 DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
334 $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
335 $ SALFAI, SALFAR, SBETA, SCALE, TEMP
336* ..
337* .. Local Arrays ..
338 LOGICAL LDUMMA( 1 )
339* ..
340* .. External Subroutines ..
341 EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
342 $ dlascl, dlaset, dorgqr, dormqr, dtgevc, xerbla
343* ..
344* .. External Functions ..
345 LOGICAL LSAME
346 INTEGER ILAENV
347 DOUBLE PRECISION DLAMCH, DLANGE
348 EXTERNAL lsame, ilaenv, dlamch, dlange
349* ..
350* .. Intrinsic Functions ..
351 INTRINSIC abs, int, max
352* ..
353* .. Executable Statements ..
354*
355* Decode the input arguments
356*
357 IF( lsame( jobvl, 'N' ) ) THEN
358 ijobvl = 1
359 ilvl = .false.
360 ELSE IF( lsame( jobvl, 'V' ) ) THEN
361 ijobvl = 2
362 ilvl = .true.
363 ELSE
364 ijobvl = -1
365 ilvl = .false.
366 END IF
367*
368 IF( lsame( jobvr, 'N' ) ) THEN
369 ijobvr = 1
370 ilvr = .false.
371 ELSE IF( lsame( jobvr, 'V' ) ) THEN
372 ijobvr = 2
373 ilvr = .true.
374 ELSE
375 ijobvr = -1
376 ilvr = .false.
377 END IF
378 ilv = ilvl .OR. ilvr
379*
380* Test the input arguments
381*
382 lwkmin = max( 8*n, 1 )
383 lwkopt = lwkmin
384 work( 1 ) = lwkopt
385 lquery = ( lwork.EQ.-1 )
386 info = 0
387 IF( ijobvl.LE.0 ) THEN
388 info = -1
389 ELSE IF( ijobvr.LE.0 ) THEN
390 info = -2
391 ELSE IF( n.LT.0 ) THEN
392 info = -3
393 ELSE IF( lda.LT.max( 1, n ) ) THEN
394 info = -5
395 ELSE IF( ldb.LT.max( 1, n ) ) THEN
396 info = -7
397 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
398 info = -12
399 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
400 info = -14
401 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
402 info = -16
403 END IF
404*
405 IF( info.EQ.0 ) THEN
406 nb1 = ilaenv( 1, 'DGEQRF', ' ', n, n, -1, -1 )
407 nb2 = ilaenv( 1, 'DORMQR', ' ', n, n, n, -1 )
408 nb3 = ilaenv( 1, 'DORGQR', ' ', n, n, n, -1 )
409 nb = max( nb1, nb2, nb3 )
410 lopt = 2*n + max( 6*n, n*( nb+1 ) )
411 work( 1 ) = lopt
412 END IF
413*
414 IF( info.NE.0 ) THEN
415 CALL xerbla( 'DGEGV ', -info )
416 RETURN
417 ELSE IF( lquery ) THEN
418 RETURN
419 END IF
420*
421* Quick return if possible
422*
423 IF( n.EQ.0 )
424 $ RETURN
425*
426* Get machine constants
427*
428 eps = dlamch( 'E' )*dlamch( 'B' )
429 safmin = dlamch( 'S' )
430 safmin = safmin + safmin
431 safmax = one / safmin
432 onepls = one + ( 4*eps )
433*
434* Scale A
435*
436 anrm = dlange( 'M', n, n, a, lda, work )
437 anrm1 = anrm
438 anrm2 = one
439 IF( anrm.LT.one ) THEN
440 IF( safmax*anrm.LT.one ) THEN
441 anrm1 = safmin
442 anrm2 = safmax*anrm
443 END IF
444 END IF
445*
446 IF( anrm.GT.zero ) THEN
447 CALL dlascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
448 IF( iinfo.NE.0 ) THEN
449 info = n + 10
450 RETURN
451 END IF
452 END IF
453*
454* Scale B
455*
456 bnrm = dlange( 'M', n, n, b, ldb, work )
457 bnrm1 = bnrm
458 bnrm2 = one
459 IF( bnrm.LT.one ) THEN
460 IF( safmax*bnrm.LT.one ) THEN
461 bnrm1 = safmin
462 bnrm2 = safmax*bnrm
463 END IF
464 END IF
465*
466 IF( bnrm.GT.zero ) THEN
467 CALL dlascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
468 IF( iinfo.NE.0 ) THEN
469 info = n + 10
470 RETURN
471 END IF
472 END IF
473*
474* Permute the matrix to make it more nearly triangular
475* Workspace layout: (8*N words -- "work" requires 6*N words)
476* left_permutation, right_permutation, work...
477*
478 ileft = 1
479 iright = n + 1
480 iwork = iright + n
481 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
482 $ work( iright ), work( iwork ), iinfo )
483 IF( iinfo.NE.0 ) THEN
484 info = n + 1
485 GO TO 120
486 END IF
487*
488* Reduce B to triangular form, and initialize VL and/or VR
489* Workspace layout: ("work..." must have at least N words)
490* left_permutation, right_permutation, tau, work...
491*
492 irows = ihi + 1 - ilo
493 IF( ilv ) THEN
494 icols = n + 1 - ilo
495 ELSE
496 icols = irows
497 END IF
498 itau = iwork
499 iwork = itau + irows
500 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
501 $ work( iwork ), lwork+1-iwork, iinfo )
502 IF( iinfo.GE.0 )
503 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
504 IF( iinfo.NE.0 ) THEN
505 info = n + 2
506 GO TO 120
507 END IF
508*
509 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
510 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
511 $ lwork+1-iwork, iinfo )
512 IF( iinfo.GE.0 )
513 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
514 IF( iinfo.NE.0 ) THEN
515 info = n + 3
516 GO TO 120
517 END IF
518*
519 IF( ilvl ) THEN
520 CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
521 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
522 $ vl( ilo+1, ilo ), ldvl )
523 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
524 $ work( itau ), work( iwork ), lwork+1-iwork,
525 $ iinfo )
526 IF( iinfo.GE.0 )
527 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
528 IF( iinfo.NE.0 ) THEN
529 info = n + 4
530 GO TO 120
531 END IF
532 END IF
533*
534 IF( ilvr )
535 $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
536*
537* Reduce to generalized Hessenberg form
538*
539 IF( ilv ) THEN
540*
541* Eigenvectors requested -- work on whole matrix.
542*
543 CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
544 $ ldvl, vr, ldvr, iinfo )
545 ELSE
546 CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
547 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
548 END IF
549 IF( iinfo.NE.0 ) THEN
550 info = n + 5
551 GO TO 120
552 END IF
553*
554* Perform QZ algorithm
555* Workspace layout: ("work..." must have at least 1 word)
556* left_permutation, right_permutation, work...
557*
558 iwork = itau
559 IF( ilv ) THEN
560 chtemp = 'S'
561 ELSE
562 chtemp = 'E'
563 END IF
564 CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
565 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
566 $ work( iwork ), lwork+1-iwork, iinfo )
567 IF( iinfo.GE.0 )
568 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
569 IF( iinfo.NE.0 ) THEN
570 IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
571 info = iinfo
572 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
573 info = iinfo - n
574 ELSE
575 info = n + 6
576 END IF
577 GO TO 120
578 END IF
579*
580 IF( ilv ) THEN
581*
582* Compute Eigenvectors (DTGEVC requires 6*N words of workspace)
583*
584 IF( ilvl ) THEN
585 IF( ilvr ) THEN
586 chtemp = 'B'
587 ELSE
588 chtemp = 'L'
589 END IF
590 ELSE
591 chtemp = 'R'
592 END IF
593*
594 CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
595 $ vr, ldvr, n, in, work( iwork ), iinfo )
596 IF( iinfo.NE.0 ) THEN
597 info = n + 7
598 GO TO 120
599 END IF
600*
601* Undo balancing on VL and VR, rescale
602*
603 IF( ilvl ) THEN
604 CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
605 $ work( iright ), n, vl, ldvl, iinfo )
606 IF( iinfo.NE.0 ) THEN
607 info = n + 8
608 GO TO 120
609 END IF
610 DO 50 jc = 1, n
611 IF( alphai( jc ).LT.zero )
612 $ GO TO 50
613 temp = zero
614 IF( alphai( jc ).EQ.zero ) THEN
615 DO 10 jr = 1, n
616 temp = max( temp, abs( vl( jr, jc ) ) )
617 10 CONTINUE
618 ELSE
619 DO 20 jr = 1, n
620 temp = max( temp, abs( vl( jr, jc ) )+
621 $ abs( vl( jr, jc+1 ) ) )
622 20 CONTINUE
623 END IF
624 IF( temp.LT.safmin )
625 $ GO TO 50
626 temp = one / temp
627 IF( alphai( jc ).EQ.zero ) THEN
628 DO 30 jr = 1, n
629 vl( jr, jc ) = vl( jr, jc )*temp
630 30 CONTINUE
631 ELSE
632 DO 40 jr = 1, n
633 vl( jr, jc ) = vl( jr, jc )*temp
634 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
635 40 CONTINUE
636 END IF
637 50 CONTINUE
638 END IF
639 IF( ilvr ) THEN
640 CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
641 $ work( iright ), n, vr, ldvr, iinfo )
642 IF( iinfo.NE.0 ) THEN
643 info = n + 9
644 GO TO 120
645 END IF
646 DO 100 jc = 1, n
647 IF( alphai( jc ).LT.zero )
648 $ GO TO 100
649 temp = zero
650 IF( alphai( jc ).EQ.zero ) THEN
651 DO 60 jr = 1, n
652 temp = max( temp, abs( vr( jr, jc ) ) )
653 60 CONTINUE
654 ELSE
655 DO 70 jr = 1, n
656 temp = max( temp, abs( vr( jr, jc ) )+
657 $ abs( vr( jr, jc+1 ) ) )
658 70 CONTINUE
659 END IF
660 IF( temp.LT.safmin )
661 $ GO TO 100
662 temp = one / temp
663 IF( alphai( jc ).EQ.zero ) THEN
664 DO 80 jr = 1, n
665 vr( jr, jc ) = vr( jr, jc )*temp
666 80 CONTINUE
667 ELSE
668 DO 90 jr = 1, n
669 vr( jr, jc ) = vr( jr, jc )*temp
670 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
671 90 CONTINUE
672 END IF
673 100 CONTINUE
674 END IF
675*
676* End of eigenvector calculation
677*
678 END IF
679*
680* Undo scaling in alpha, beta
681*
682* Note: this does not give the alpha and beta for the unscaled
683* problem.
684*
685* Un-scaling is limited to avoid underflow in alpha and beta
686* if they are significant.
687*
688 DO 110 jc = 1, n
689 absar = abs( alphar( jc ) )
690 absai = abs( alphai( jc ) )
691 absb = abs( beta( jc ) )
692 salfar = anrm*alphar( jc )
693 salfai = anrm*alphai( jc )
694 sbeta = bnrm*beta( jc )
695 ilimit = .false.
696 scale = one
697*
698* Check for significant underflow in ALPHAI
699*
700 IF( abs( salfai ).LT.safmin .AND. absai.GE.
701 $ max( safmin, eps*absar, eps*absb ) ) THEN
702 ilimit = .true.
703 scale = ( onepls*safmin / anrm1 ) /
704 $ max( onepls*safmin, anrm2*absai )
705*
706 ELSE IF( salfai.EQ.zero ) THEN
707*
708* If insignificant underflow in ALPHAI, then make the
709* conjugate eigenvalue real.
710*
711 IF( alphai( jc ).LT.zero .AND. jc.GT.1 ) THEN
712 alphai( jc-1 ) = zero
713 ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n ) THEN
714 alphai( jc+1 ) = zero
715 END IF
716 END IF
717*
718* Check for significant underflow in ALPHAR
719*
720 IF( abs( salfar ).LT.safmin .AND. absar.GE.
721 $ max( safmin, eps*absai, eps*absb ) ) THEN
722 ilimit = .true.
723 scale = max( scale, ( onepls*safmin / anrm1 ) /
724 $ max( onepls*safmin, anrm2*absar ) )
725 END IF
726*
727* Check for significant underflow in BETA
728*
729 IF( abs( sbeta ).LT.safmin .AND. absb.GE.
730 $ max( safmin, eps*absar, eps*absai ) ) THEN
731 ilimit = .true.
732 scale = max( scale, ( onepls*safmin / bnrm1 ) /
733 $ max( onepls*safmin, bnrm2*absb ) )
734 END IF
735*
736* Check for possible overflow when limiting scaling
737*
738 IF( ilimit ) THEN
739 temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
740 $ abs( sbeta ) )
741 IF( temp.GT.one )
742 $ scale = scale / temp
743 IF( scale.LT.one )
744 $ ilimit = .false.
745 END IF
746*
747* Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
748*
749 IF( ilimit ) THEN
750 salfar = ( scale*alphar( jc ) )*anrm
751 salfai = ( scale*alphai( jc ) )*anrm
752 sbeta = ( scale*beta( jc ) )*bnrm
753 END IF
754 alphar( jc ) = salfar
755 alphai( jc ) = salfai
756 beta( jc ) = sbeta
757 110 CONTINUE
758*
759 120 CONTINUE
760 work( 1 ) = lwkopt
761*
762 RETURN
763*
764* End of DGEGV
765*
dtgevc
subroutine dtgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
DTGEVC
Definition dtgevc.f:295
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339

◆ dgges()

subroutine dgges ( character jobvsl,
character jobvsr,
character sort,
external selctg,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
integer sdim,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvsl, * ) vsl,
integer ldvsl,
double precision, dimension( ldvsr, * ) vsr,
integer ldvsr,
double precision, dimension( * ) work,
integer lwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
!> the generalized eigenvalues, the generalized real Schur form (S,T),
!> optionally, the left and/or right matrices of Schur vectors (VSL and
!> VSR). This gives the generalized Schur factorization
!>
!>          (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-triangular matrix S and the upper triangular matrix T.The
!> leading columns of VSL and VSR then form an orthonormal basis for the
!> corresponding left and right eigenspaces (deflating subspaces).
!>
!> (If only the generalized eigenvalues are needed, use the driver
!> DGGEV 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 or both being zero.
!>
!> A pair of matrices (S,T) is in generalized real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!>  by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!>
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 three DOUBLE PRECISION 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 (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>
!>          Note that in the ill-conditioned case, a selected complex
!>          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
!>          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
!>          in this case.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION 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 DOUBLE PRECISION 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.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
!>          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N = 0, LWORK >= 1, else LWORK >= 8*N+16.
!>          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]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 ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
!>                =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 DTGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 281 of file dgges.f.

284*
285* -- LAPACK driver routine --
286* -- LAPACK is a software package provided by Univ. of Tennessee, --
287* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288*
289* .. Scalar Arguments ..
290 CHARACTER JOBVSL, JOBVSR, SORT
291 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
292* ..
293* .. Array Arguments ..
294 LOGICAL BWORK( * )
295 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
296 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
297 $ VSR( LDVSR, * ), WORK( * )
298* ..
299* .. Function Arguments ..
300 LOGICAL SELCTG
301 EXTERNAL selctg
302* ..
303*
304* =====================================================================
305*
306* .. Parameters ..
307 DOUBLE PRECISION ZERO, ONE
308 parameter( zero = 0.0d+0, one = 1.0d+0 )
309* ..
310* .. Local Scalars ..
311 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
312 $ LQUERY, LST2SL, WANTST
313 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
314 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
315 $ MINWRK
316 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
317 $ PVSR, SAFMAX, SAFMIN, SMLNUM
318* ..
319* .. Local Arrays ..
320 INTEGER IDUM( 1 )
321 DOUBLE PRECISION DIF( 2 )
322* ..
323* .. External Subroutines ..
324 EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
325 $ dlacpy, dlascl, dlaset, dorgqr, dormqr, dtgsen,
326 $ xerbla
327* ..
328* .. External Functions ..
329 LOGICAL LSAME
330 INTEGER ILAENV
331 DOUBLE PRECISION DLAMCH, DLANGE
332 EXTERNAL lsame, ilaenv, dlamch, dlange
333* ..
334* .. Intrinsic Functions ..
335 INTRINSIC abs, max, sqrt
336* ..
337* .. Executable Statements ..
338*
339* Decode the input arguments
340*
341 IF( lsame( jobvsl, 'N' ) ) THEN
342 ijobvl = 1
343 ilvsl = .false.
344 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
345 ijobvl = 2
346 ilvsl = .true.
347 ELSE
348 ijobvl = -1
349 ilvsl = .false.
350 END IF
351*
352 IF( lsame( jobvsr, 'N' ) ) THEN
353 ijobvr = 1
354 ilvsr = .false.
355 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
356 ijobvr = 2
357 ilvsr = .true.
358 ELSE
359 ijobvr = -1
360 ilvsr = .false.
361 END IF
362*
363 wantst = lsame( sort, 'S' )
364*
365* Test the input arguments
366*
367 info = 0
368 lquery = ( lwork.EQ.-1 )
369 IF( ijobvl.LE.0 ) THEN
370 info = -1
371 ELSE IF( ijobvr.LE.0 ) THEN
372 info = -2
373 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
374 info = -3
375 ELSE IF( n.LT.0 ) THEN
376 info = -5
377 ELSE IF( lda.LT.max( 1, n ) ) THEN
378 info = -7
379 ELSE IF( ldb.LT.max( 1, n ) ) THEN
380 info = -9
381 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
382 info = -15
383 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
384 info = -17
385 END IF
386*
387* Compute workspace
388* (Note: Comments in the code beginning "Workspace:" describe the
389* minimal amount of workspace needed at that point in the code,
390* as well as the preferred amount for good performance.
391* NB refers to the optimal block size for the immediately
392* following subroutine, as returned by ILAENV.)
393*
394 IF( info.EQ.0 ) THEN
395 IF( n.GT.0 )THEN
396 minwrk = max( 8*n, 6*n + 16 )
397 maxwrk = minwrk - n +
398 $ n*ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 )
399 maxwrk = max( maxwrk, minwrk - n +
400 $ n*ilaenv( 1, 'DORMQR', ' ', n, 1, n, -1 ) )
401 IF( ilvsl ) THEN
402 maxwrk = max( maxwrk, minwrk - n +
403 $ n*ilaenv( 1, 'DORGQR', ' ', n, 1, n, -1 ) )
404 END IF
405 ELSE
406 minwrk = 1
407 maxwrk = 1
408 END IF
409 work( 1 ) = maxwrk
410*
411 IF( lwork.LT.minwrk .AND. .NOT.lquery )
412 $ info = -19
413 END IF
414*
415 IF( info.NE.0 ) THEN
416 CALL xerbla( 'DGGES ', -info )
417 RETURN
418 ELSE IF( lquery ) THEN
419 RETURN
420 END IF
421*
422* Quick return if possible
423*
424 IF( n.EQ.0 ) THEN
425 sdim = 0
426 RETURN
427 END IF
428*
429* Get machine constants
430*
431 eps = dlamch( 'P' )
432 safmin = dlamch( 'S' )
433 safmax = one / safmin
434 CALL dlabad( safmin, safmax )
435 smlnum = sqrt( safmin ) / eps
436 bignum = one / smlnum
437*
438* Scale A if max element outside range [SMLNUM,BIGNUM]
439*
440 anrm = dlange( 'M', n, n, a, lda, work )
441 ilascl = .false.
442 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
443 anrmto = smlnum
444 ilascl = .true.
445 ELSE IF( anrm.GT.bignum ) THEN
446 anrmto = bignum
447 ilascl = .true.
448 END IF
449 IF( ilascl )
450 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
451*
452* Scale B if max element outside range [SMLNUM,BIGNUM]
453*
454 bnrm = dlange( 'M', n, n, b, ldb, work )
455 ilbscl = .false.
456 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
457 bnrmto = smlnum
458 ilbscl = .true.
459 ELSE IF( bnrm.GT.bignum ) THEN
460 bnrmto = bignum
461 ilbscl = .true.
462 END IF
463 IF( ilbscl )
464 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
465*
466* Permute the matrix to make it more nearly triangular
467* (Workspace: need 6*N + 2*N space for storing balancing factors)
468*
469 ileft = 1
470 iright = n + 1
471 iwrk = iright + n
472 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
473 $ work( iright ), work( iwrk ), ierr )
474*
475* Reduce B to triangular form (QR decomposition of B)
476* (Workspace: need N, prefer N*NB)
477*
478 irows = ihi + 1 - ilo
479 icols = n + 1 - ilo
480 itau = iwrk
481 iwrk = itau + irows
482 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
483 $ work( iwrk ), lwork+1-iwrk, ierr )
484*
485* Apply the orthogonal transformation to matrix A
486* (Workspace: need N, prefer N*NB)
487*
488 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
489 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
490 $ lwork+1-iwrk, ierr )
491*
492* Initialize VSL
493* (Workspace: need N, prefer N*NB)
494*
495 IF( ilvsl ) THEN
496 CALL dlaset( 'Full', n, n, zero, one, vsl, ldvsl )
497 IF( irows.GT.1 ) THEN
498 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
499 $ vsl( ilo+1, ilo ), ldvsl )
500 END IF
501 CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
502 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
503 END IF
504*
505* Initialize VSR
506*
507 IF( ilvsr )
508 $ CALL dlaset( 'Full', n, n, zero, one, vsr, ldvsr )
509*
510* Reduce to generalized Hessenberg form
511* (Workspace: none needed)
512*
513 CALL dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
514 $ ldvsl, vsr, ldvsr, ierr )
515*
516* Perform QZ algorithm, computing Schur vectors if desired
517* (Workspace: need N)
518*
519 iwrk = itau
520 CALL dhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
521 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
522 $ work( iwrk ), lwork+1-iwrk, ierr )
523 IF( ierr.NE.0 ) THEN
524 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
525 info = ierr
526 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
527 info = ierr - n
528 ELSE
529 info = n + 1
530 END IF
531 GO TO 50
532 END IF
533*
534* Sort eigenvalues ALPHA/BETA if desired
535* (Workspace: need 4*N+16 )
536*
537 sdim = 0
538 IF( wantst ) THEN
539*
540* Undo scaling on eigenvalues before SELCTGing
541*
542 IF( ilascl ) THEN
543 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
544 $ ierr )
545 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
546 $ ierr )
547 END IF
548 IF( ilbscl )
549 $ CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
550*
551* Select eigenvalues
552*
553 DO 10 i = 1, n
554 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
555 10 CONTINUE
556*
557 CALL dtgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
558 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
559 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
560 $ ierr )
561 IF( ierr.EQ.1 )
562 $ info = n + 3
563*
564 END IF
565*
566* Apply back-permutation to VSL and VSR
567* (Workspace: none needed)
568*
569 IF( ilvsl )
570 $ CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
571 $ work( iright ), n, vsl, ldvsl, ierr )
572*
573 IF( ilvsr )
574 $ CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
575 $ work( iright ), n, vsr, ldvsr, ierr )
576*
577* Check if unscaling would cause over/underflow, if so, rescale
578* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
579* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
580*
581 IF( ilascl ) THEN
582 DO 20 i = 1, n
583 IF( alphai( i ).NE.zero ) THEN
584 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
585 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) ) THEN
586 work( 1 ) = abs( a( i, i ) / alphar( i ) )
587 beta( i ) = beta( i )*work( 1 )
588 alphar( i ) = alphar( i )*work( 1 )
589 alphai( i ) = alphai( i )*work( 1 )
590 ELSE IF( ( alphai( i ) / safmax ).GT.
591 $ ( anrmto / anrm ) .OR.
592 $ ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
593 $ THEN
594 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
595 beta( i ) = beta( i )*work( 1 )
596 alphar( i ) = alphar( i )*work( 1 )
597 alphai( i ) = alphai( i )*work( 1 )
598 END IF
599 END IF
600 20 CONTINUE
601 END IF
602*
603 IF( ilbscl ) THEN
604 DO 30 i = 1, n
605 IF( alphai( i ).NE.zero ) THEN
606 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
607 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) ) THEN
608 work( 1 ) = abs( b( i, i ) / beta( i ) )
609 beta( i ) = beta( i )*work( 1 )
610 alphar( i ) = alphar( i )*work( 1 )
611 alphai( i ) = alphai( i )*work( 1 )
612 END IF
613 END IF
614 30 CONTINUE
615 END IF
616*
617* Undo scaling
618*
619 IF( ilascl ) THEN
620 CALL dlascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
621 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
622 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
623 END IF
624*
625 IF( ilbscl ) THEN
626 CALL dlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
627 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
628 END IF
629*
630 IF( wantst ) THEN
631*
632* Check if reordering is correct
633*
634 lastsl = .true.
635 lst2sl = .true.
636 sdim = 0
637 ip = 0
638 DO 40 i = 1, n
639 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
640 IF( alphai( i ).EQ.zero ) THEN
641 IF( cursl )
642 $ sdim = sdim + 1
643 ip = 0
644 IF( cursl .AND. .NOT.lastsl )
645 $ info = n + 2
646 ELSE
647 IF( ip.EQ.1 ) THEN
648*
649* Last eigenvalue of conjugate pair
650*
651 cursl = cursl .OR. lastsl
652 lastsl = cursl
653 IF( cursl )
654 $ sdim = sdim + 2
655 ip = -1
656 IF( cursl .AND. .NOT.lst2sl )
657 $ info = n + 2
658 ELSE
659*
660* First eigenvalue of conjugate pair
661*
662 ip = 1
663 END IF
664 END IF
665 lst2sl = lastsl
666 lastsl = cursl
667 40 CONTINUE
668*
669 END IF
670*
671 50 CONTINUE
672*
673 work( 1 ) = maxwrk
674*
675 RETURN
676*
677* End of DGGES
678*
dtgsen
subroutine dtgsen(ijob, wantq, wantz, select, n, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, m, pl, pr, dif, work, lwork, iwork, liwork, info)
DTGSEN
Definition dtgsen.f:451

◆ dgges3()

subroutine dgges3 ( character jobvsl,
character jobvsr,
character sort,
external selctg,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
integer sdim,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvsl, * ) vsl,
integer ldvsl,
double precision, dimension( ldvsr, * ) vsr,
integer ldvsr,
double precision, dimension( * ) work,
integer lwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> DGGES3 computes for a pair of N-by-N real nonsymmetric matrices (A,B),
!> the generalized eigenvalues, the generalized real Schur form (S,T),
!> optionally, the left and/or right matrices of Schur vectors (VSL and
!> VSR). This gives the generalized Schur factorization
!>
!>          (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-triangular matrix S and the upper triangular matrix T.The
!> leading columns of VSL and VSR then form an orthonormal basis for the
!> corresponding left and right eigenspaces (deflating subspaces).
!>
!> (If only the generalized eigenvalues are needed, use the driver
!> DGGEV 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 or both being zero.
!>
!> A pair of matrices (S,T) is in generalized real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!>  by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!>
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 three DOUBLE PRECISION 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 (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>
!>          Note that in the ill-conditioned case, a selected complex
!>          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
!>          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
!>          in this case.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION 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 DOUBLE PRECISION 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.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
!>          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>
!>          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]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 ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DLAQZ0.
!>                =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 DTGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 279 of file dgges3.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 JOBVSL, JOBVSR, SORT
289 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
290* ..
291* .. Array Arguments ..
292 LOGICAL BWORK( * )
293 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
294 $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
295 $ VSR( LDVSR, * ), WORK( * )
296* ..
297* .. Function Arguments ..
298 LOGICAL SELCTG
299 EXTERNAL selctg
300* ..
301*
302* =====================================================================
303*
304* .. Parameters ..
305 DOUBLE PRECISION ZERO, ONE
306 parameter( zero = 0.0d+0, one = 1.0d+0 )
307* ..
308* .. Local Scalars ..
309 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
310 $ LQUERY, LST2SL, WANTST
311 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
312 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, LWKOPT
313 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
314 $ PVSR, SAFMAX, SAFMIN, SMLNUM
315* ..
316* .. Local Arrays ..
317 INTEGER IDUM( 1 )
318 DOUBLE PRECISION DIF( 2 )
319* ..
320* .. External Subroutines ..
321 EXTERNAL dgeqrf, dggbak, dggbal, dgghd3, dlaqz0, dlabad,
322 $ dlacpy, dlascl, dlaset, dorgqr, dormqr, dtgsen,
323 $ xerbla
324* ..
325* .. External Functions ..
326 LOGICAL LSAME
327 DOUBLE PRECISION DLAMCH, DLANGE
328 EXTERNAL lsame, dlamch, dlange
329* ..
330* .. Intrinsic Functions ..
331 INTRINSIC abs, max, sqrt
332* ..
333* .. Executable Statements ..
334*
335* Decode the input arguments
336*
337 IF( lsame( jobvsl, 'N' ) ) THEN
338 ijobvl = 1
339 ilvsl = .false.
340 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
341 ijobvl = 2
342 ilvsl = .true.
343 ELSE
344 ijobvl = -1
345 ilvsl = .false.
346 END IF
347*
348 IF( lsame( jobvsr, 'N' ) ) THEN
349 ijobvr = 1
350 ilvsr = .false.
351 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
352 ijobvr = 2
353 ilvsr = .true.
354 ELSE
355 ijobvr = -1
356 ilvsr = .false.
357 END IF
358*
359 wantst = lsame( sort, 'S' )
360*
361* Test the input arguments
362*
363 info = 0
364 lquery = ( lwork.EQ.-1 )
365 IF( ijobvl.LE.0 ) THEN
366 info = -1
367 ELSE IF( ijobvr.LE.0 ) THEN
368 info = -2
369 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
370 info = -3
371 ELSE IF( n.LT.0 ) THEN
372 info = -5
373 ELSE IF( lda.LT.max( 1, n ) ) THEN
374 info = -7
375 ELSE IF( ldb.LT.max( 1, n ) ) THEN
376 info = -9
377 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
378 info = -15
379 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
380 info = -17
381 ELSE IF( lwork.LT.6*n+16 .AND. .NOT.lquery ) THEN
382 info = -19
383 END IF
384*
385* Compute workspace
386*
387 IF( info.EQ.0 ) THEN
388 CALL dgeqrf( n, n, b, ldb, work, work, -1, ierr )
389 lwkopt = max( 6*n+16, 3*n+int( work( 1 ) ) )
390 CALL dormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work,
391 $ -1, ierr )
392 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
393 IF( ilvsl ) THEN
394 CALL dorgqr( n, n, n, vsl, ldvsl, work, work, -1, ierr )
395 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
396 END IF
397 CALL dgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
398 $ ldvsl, vsr, ldvsr, work, -1, ierr )
399 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
400 CALL dlaqz0( 'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
401 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
402 $ work, -1, 0, ierr )
403 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
404 IF( wantst ) THEN
405 CALL dtgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
406 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
407 $ sdim, pvsl, pvsr, dif, work, -1, idum, 1,
408 $ ierr )
409 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
410 END IF
411 work( 1 ) = lwkopt
412 END IF
413*
414 IF( info.NE.0 ) THEN
415 CALL xerbla( 'DGGES3 ', -info )
416 RETURN
417 ELSE IF( lquery ) THEN
418 RETURN
419 END IF
420*
421* Quick return if possible
422*
423 IF( n.EQ.0 ) THEN
424 sdim = 0
425 RETURN
426 END IF
427*
428* Get machine constants
429*
430 eps = dlamch( 'P' )
431 safmin = dlamch( 'S' )
432 safmax = one / safmin
433 CALL dlabad( safmin, safmax )
434 smlnum = sqrt( safmin ) / eps
435 bignum = one / smlnum
436*
437* Scale A if max element outside range [SMLNUM,BIGNUM]
438*
439 anrm = dlange( 'M', n, n, a, lda, work )
440 ilascl = .false.
441 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
442 anrmto = smlnum
443 ilascl = .true.
444 ELSE IF( anrm.GT.bignum ) THEN
445 anrmto = bignum
446 ilascl = .true.
447 END IF
448 IF( ilascl )
449 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
450*
451* Scale B if max element outside range [SMLNUM,BIGNUM]
452*
453 bnrm = dlange( 'M', n, n, b, ldb, work )
454 ilbscl = .false.
455 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
456 bnrmto = smlnum
457 ilbscl = .true.
458 ELSE IF( bnrm.GT.bignum ) THEN
459 bnrmto = bignum
460 ilbscl = .true.
461 END IF
462 IF( ilbscl )
463 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
464*
465* Permute the matrix to make it more nearly triangular
466*
467 ileft = 1
468 iright = n + 1
469 iwrk = iright + n
470 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
471 $ work( iright ), work( iwrk ), ierr )
472*
473* Reduce B to triangular form (QR decomposition of B)
474*
475 irows = ihi + 1 - ilo
476 icols = n + 1 - ilo
477 itau = iwrk
478 iwrk = itau + irows
479 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
480 $ work( iwrk ), lwork+1-iwrk, ierr )
481*
482* Apply the orthogonal transformation to matrix A
483*
484 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
485 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
486 $ lwork+1-iwrk, ierr )
487*
488* Initialize VSL
489*
490 IF( ilvsl ) THEN
491 CALL dlaset( 'Full', n, n, zero, one, vsl, ldvsl )
492 IF( irows.GT.1 ) THEN
493 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
494 $ vsl( ilo+1, ilo ), ldvsl )
495 END IF
496 CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
497 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
498 END IF
499*
500* Initialize VSR
501*
502 IF( ilvsr )
503 $ CALL dlaset( 'Full', n, n, zero, one, vsr, ldvsr )
504*
505* Reduce to generalized Hessenberg form
506*
507 CALL dgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
508 $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk,
509 $ ierr )
510*
511* Perform QZ algorithm, computing Schur vectors if desired
512*
513 iwrk = itau
514 CALL dlaqz0( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
515 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
516 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
517 IF( ierr.NE.0 ) THEN
518 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
519 info = ierr
520 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
521 info = ierr - n
522 ELSE
523 info = n + 1
524 END IF
525 GO TO 50
526 END IF
527*
528* Sort eigenvalues ALPHA/BETA if desired
529*
530 sdim = 0
531 IF( wantst ) THEN
532*
533* Undo scaling on eigenvalues before SELCTGing
534*
535 IF( ilascl ) THEN
536 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
537 $ ierr )
538 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
539 $ ierr )
540 END IF
541 IF( ilbscl )
542 $ CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
543*
544* Select eigenvalues
545*
546 DO 10 i = 1, n
547 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
548 10 CONTINUE
549*
550 CALL dtgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
551 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
552 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
553 $ ierr )
554 IF( ierr.EQ.1 )
555 $ info = n + 3
556*
557 END IF
558*
559* Apply back-permutation to VSL and VSR
560*
561 IF( ilvsl )
562 $ CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
563 $ work( iright ), n, vsl, ldvsl, ierr )
564*
565 IF( ilvsr )
566 $ CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
567 $ work( iright ), n, vsr, ldvsr, ierr )
568*
569* Check if unscaling would cause over/underflow, if so, rescale
570* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
571* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
572*
573 IF( ilascl ) THEN
574 DO 20 i = 1, n
575 IF( alphai( i ).NE.zero ) THEN
576 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
577 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) ) THEN
578 work( 1 ) = abs( a( i, i ) / alphar( i ) )
579 beta( i ) = beta( i )*work( 1 )
580 alphar( i ) = alphar( i )*work( 1 )
581 alphai( i ) = alphai( i )*work( 1 )
582 ELSE IF( ( alphai( i ) / safmax ).GT.
583 $ ( anrmto / anrm ) .OR.
584 $ ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
585 $ THEN
586 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
587 beta( i ) = beta( i )*work( 1 )
588 alphar( i ) = alphar( i )*work( 1 )
589 alphai( i ) = alphai( i )*work( 1 )
590 END IF
591 END IF
592 20 CONTINUE
593 END IF
594*
595 IF( ilbscl ) THEN
596 DO 30 i = 1, n
597 IF( alphai( i ).NE.zero ) THEN
598 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
599 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) ) THEN
600 work( 1 ) = abs( b( i, i ) / beta( i ) )
601 beta( i ) = beta( i )*work( 1 )
602 alphar( i ) = alphar( i )*work( 1 )
603 alphai( i ) = alphai( i )*work( 1 )
604 END IF
605 END IF
606 30 CONTINUE
607 END IF
608*
609* Undo scaling
610*
611 IF( ilascl ) THEN
612 CALL dlascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
613 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
614 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
615 END IF
616*
617 IF( ilbscl ) THEN
618 CALL dlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
619 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
620 END IF
621*
622 IF( wantst ) THEN
623*
624* Check if reordering is correct
625*
626 lastsl = .true.
627 lst2sl = .true.
628 sdim = 0
629 ip = 0
630 DO 40 i = 1, n
631 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
632 IF( alphai( i ).EQ.zero ) THEN
633 IF( cursl )
634 $ sdim = sdim + 1
635 ip = 0
636 IF( cursl .AND. .NOT.lastsl )
637 $ info = n + 2
638 ELSE
639 IF( ip.EQ.1 ) THEN
640*
641* Last eigenvalue of conjugate pair
642*
643 cursl = cursl .OR. lastsl
644 lastsl = cursl
645 IF( cursl )
646 $ sdim = sdim + 2
647 ip = -1
648 IF( cursl .AND. .NOT.lst2sl )
649 $ info = n + 2
650 ELSE
651*
652* First eigenvalue of conjugate pair
653*
654 ip = 1
655 END IF
656 END IF
657 lst2sl = lastsl
658 lastsl = cursl
659 40 CONTINUE
660*
661 END IF
662*
663 50 CONTINUE
664*
665 work( 1 ) = lwkopt
666*
667 RETURN
668*
669* End of DGGES3
670*
dlaqz0
recursive subroutine dlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, rec, info)
DLAQZ0
Definition dlaqz0.f:306
dgghd3
subroutine dgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
DGGHD3
Definition dgghd3.f:230

◆ dggesx()

subroutine dggesx ( character jobvsl,
character jobvsr,
character sort,
external selctg,
character sense,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
integer sdim,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvsl, * ) vsl,
integer ldvsl,
double precision, dimension( ldvsr, * ) vsr,
integer ldvsr,
double precision, dimension( 2 ) rconde,
double precision, dimension( 2 ) rcondv,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> DGGESX computes for a pair of N-by-N real nonsymmetric matrices
!> (A,B), the generalized eigenvalues, the real 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)**T, (VSL) T (VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-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 real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!>  by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!>
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 three DOUBLE PRECISION 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 (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>          Note that a selected complex eigenvalue may no longer satisfy
!>          SELCTG(ALPHAR(j),ALPHAI(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.
!>
[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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!>
[out]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i
!>          and BETA(j),j=1,...,N  are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDVSR,N)
!>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
!>          Not referenced if JOBVSR = 'N'.
!>
[in]LDVSR
!>          LDVSR is INTEGER
!>          The leading dimension of the matrix VSR. LDVSR >= 1, and
!>          if JOBVSR = 'V', LDVSR >= N.
!>
[out]RCONDE
!>          RCONDE is DOUBLE PRECISION array, dimension ( 2 )
!>          If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
!>          reciprocal condition numbers for the average of the selected
!>          eigenvalues.
!>          Not referenced if SENSE = 'N' or 'V'.
!>
[out]RCONDV
!>          RCONDV is DOUBLE PRECISION array, dimension ( 2 )
!>          If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
!>          reciprocal condition numbers for the selected deflating
!>          subspaces.
!>          Not referenced if SENSE = 'N' or 'E'.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
!>          LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
!>          LWORK >= max( 8*N, 6*N+16 ).
!>          Note that 2*SDIM*(N-SDIM) <= N*N/2.
!>          Note also that an error is only returned if
!>          LWORK < max( 8*N, 6*N+16), 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]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
!>          LIWORK >= N+6.
!>
!>          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 ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DHGEQZ
!>                =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 DTGSEN.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  An approximate (asymptotic) bound on the average absolute error of
!>  the selected eigenvalues is
!>
!>       EPS * norm((A, B)) / RCONDE( 1 ).
!>
!>  An approximate (asymptotic) bound on the maximum angular error in
!>  the computed deflating subspaces is
!>
!>       EPS * norm((A, B)) / RCONDV( 2 ).
!>
!>  See LAPACK User's Guide, section 4.11 for more information.
!>

Definition at line 361 of file dggesx.f.

365*
366* -- LAPACK driver routine --
367* -- LAPACK is a software package provided by Univ. of Tennessee, --
368* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
369*
370* .. Scalar Arguments ..
371 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
372 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
373 $ SDIM
374* ..
375* .. Array Arguments ..
376 LOGICAL BWORK( * )
377 INTEGER IWORK( * )
378 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
379 $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
380 $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
381 $ WORK( * )
382* ..
383* .. Function Arguments ..
384 LOGICAL SELCTG
385 EXTERNAL selctg
386* ..
387*
388* =====================================================================
389*
390* .. Parameters ..
391 DOUBLE PRECISION ZERO, ONE
392 parameter( zero = 0.0d+0, one = 1.0d+0 )
393* ..
394* .. Local Scalars ..
395 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
396 $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
397 $ WANTSV
398 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
399 $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
400 $ LIWMIN, LWRK, MAXWRK, MINWRK
401 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
402 $ PR, SAFMAX, SAFMIN, SMLNUM
403* ..
404* .. Local Arrays ..
405 DOUBLE PRECISION DIF( 2 )
406* ..
407* .. External Subroutines ..
408 EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
409 $ dlacpy, dlascl, dlaset, dorgqr, dormqr, dtgsen,
410 $ xerbla
411* ..
412* .. External Functions ..
413 LOGICAL LSAME
414 INTEGER ILAENV
415 DOUBLE PRECISION DLAMCH, DLANGE
416 EXTERNAL lsame, ilaenv, dlamch, dlange
417* ..
418* .. Intrinsic Functions ..
419 INTRINSIC abs, max, sqrt
420* ..
421* .. Executable Statements ..
422*
423* Decode the input arguments
424*
425 IF( lsame( jobvsl, 'N' ) ) THEN
426 ijobvl = 1
427 ilvsl = .false.
428 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
429 ijobvl = 2
430 ilvsl = .true.
431 ELSE
432 ijobvl = -1
433 ilvsl = .false.
434 END IF
435*
436 IF( lsame( jobvsr, 'N' ) ) THEN
437 ijobvr = 1
438 ilvsr = .false.
439 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
440 ijobvr = 2
441 ilvsr = .true.
442 ELSE
443 ijobvr = -1
444 ilvsr = .false.
445 END IF
446*
447 wantst = lsame( sort, 'S' )
448 wantsn = lsame( sense, 'N' )
449 wantse = lsame( sense, 'E' )
450 wantsv = lsame( sense, 'V' )
451 wantsb = lsame( sense, 'B' )
452 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
453 IF( wantsn ) THEN
454 ijob = 0
455 ELSE IF( wantse ) THEN
456 ijob = 1
457 ELSE IF( wantsv ) THEN
458 ijob = 2
459 ELSE IF( wantsb ) THEN
460 ijob = 4
461 END IF
462*
463* Test the input arguments
464*
465 info = 0
466 IF( ijobvl.LE.0 ) THEN
467 info = -1
468 ELSE IF( ijobvr.LE.0 ) THEN
469 info = -2
470 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
471 info = -3
472 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
473 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
474 info = -5
475 ELSE IF( n.LT.0 ) THEN
476 info = -6
477 ELSE IF( lda.LT.max( 1, n ) ) THEN
478 info = -8
479 ELSE IF( ldb.LT.max( 1, n ) ) THEN
480 info = -10
481 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
482 info = -16
483 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
484 info = -18
485 END IF
486*
487* Compute workspace
488* (Note: Comments in the code beginning "Workspace:" describe the
489* minimal amount of workspace needed at that point in the code,
490* as well as the preferred amount for good performance.
491* NB refers to the optimal block size for the immediately
492* following subroutine, as returned by ILAENV.)
493*
494 IF( info.EQ.0 ) THEN
495 IF( n.GT.0) THEN
496 minwrk = max( 8*n, 6*n + 16 )
497 maxwrk = minwrk - n +
498 $ n*ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 )
499 maxwrk = max( maxwrk, minwrk - n +
500 $ n*ilaenv( 1, 'DORMQR', ' ', n, 1, n, -1 ) )
501 IF( ilvsl ) THEN
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*ilaenv( 1, 'DORGQR', ' ', n, 1, n, -1 ) )
504 END IF
505 lwrk = maxwrk
506 IF( ijob.GE.1 )
507 $ lwrk = max( lwrk, n*n/2 )
508 ELSE
509 minwrk = 1
510 maxwrk = 1
511 lwrk = 1
512 END IF
513 work( 1 ) = lwrk
514 IF( wantsn .OR. n.EQ.0 ) THEN
515 liwmin = 1
516 ELSE
517 liwmin = n + 6
518 END IF
519 iwork( 1 ) = liwmin
520*
521 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
522 info = -22
523 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
524 info = -24
525 END IF
526 END IF
527*
528 IF( info.NE.0 ) THEN
529 CALL xerbla( 'DGGESX', -info )
530 RETURN
531 ELSE IF (lquery) THEN
532 RETURN
533 END IF
534*
535* Quick return if possible
536*
537 IF( n.EQ.0 ) THEN
538 sdim = 0
539 RETURN
540 END IF
541*
542* Get machine constants
543*
544 eps = dlamch( 'P' )
545 safmin = dlamch( 'S' )
546 safmax = one / safmin
547 CALL dlabad( safmin, safmax )
548 smlnum = sqrt( safmin ) / eps
549 bignum = one / smlnum
550*
551* Scale A if max element outside range [SMLNUM,BIGNUM]
552*
553 anrm = dlange( 'M', n, n, a, lda, work )
554 ilascl = .false.
555 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
556 anrmto = smlnum
557 ilascl = .true.
558 ELSE IF( anrm.GT.bignum ) THEN
559 anrmto = bignum
560 ilascl = .true.
561 END IF
562 IF( ilascl )
563 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
564*
565* Scale B if max element outside range [SMLNUM,BIGNUM]
566*
567 bnrm = dlange( 'M', n, n, b, ldb, work )
568 ilbscl = .false.
569 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
570 bnrmto = smlnum
571 ilbscl = .true.
572 ELSE IF( bnrm.GT.bignum ) THEN
573 bnrmto = bignum
574 ilbscl = .true.
575 END IF
576 IF( ilbscl )
577 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
578*
579* Permute the matrix to make it more nearly triangular
580* (Workspace: need 6*N + 2*N for permutation parameters)
581*
582 ileft = 1
583 iright = n + 1
584 iwrk = iright + n
585 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
586 $ work( iright ), work( iwrk ), ierr )
587*
588* Reduce B to triangular form (QR decomposition of B)
589* (Workspace: need N, prefer N*NB)
590*
591 irows = ihi + 1 - ilo
592 icols = n + 1 - ilo
593 itau = iwrk
594 iwrk = itau + irows
595 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
596 $ work( iwrk ), lwork+1-iwrk, ierr )
597*
598* Apply the orthogonal transformation to matrix A
599* (Workspace: need N, prefer N*NB)
600*
601 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
602 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
603 $ lwork+1-iwrk, ierr )
604*
605* Initialize VSL
606* (Workspace: need N, prefer N*NB)
607*
608 IF( ilvsl ) THEN
609 CALL dlaset( 'Full', n, n, zero, one, vsl, ldvsl )
610 IF( irows.GT.1 ) THEN
611 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
612 $ vsl( ilo+1, ilo ), ldvsl )
613 END IF
614 CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
615 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
616 END IF
617*
618* Initialize VSR
619*
620 IF( ilvsr )
621 $ CALL dlaset( 'Full', n, n, zero, one, vsr, ldvsr )
622*
623* Reduce to generalized Hessenberg form
624* (Workspace: none needed)
625*
626 CALL dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
627 $ ldvsl, vsr, ldvsr, ierr )
628*
629 sdim = 0
630*
631* Perform QZ algorithm, computing Schur vectors if desired
632* (Workspace: need N)
633*
634 iwrk = itau
635 CALL dhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
636 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
637 $ work( iwrk ), lwork+1-iwrk, ierr )
638 IF( ierr.NE.0 ) THEN
639 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
640 info = ierr
641 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
642 info = ierr - n
643 ELSE
644 info = n + 1
645 END IF
646 GO TO 60
647 END IF
648*
649* Sort eigenvalues ALPHA/BETA and compute the reciprocal of
650* condition number(s)
651* (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
652* otherwise, need 8*(N+1) )
653*
654 IF( wantst ) THEN
655*
656* Undo scaling on eigenvalues before SELCTGing
657*
658 IF( ilascl ) THEN
659 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
660 $ ierr )
661 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
662 $ ierr )
663 END IF
664 IF( ilbscl )
665 $ CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
666*
667* Select eigenvalues
668*
669 DO 10 i = 1, n
670 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
671 10 CONTINUE
672*
673* Reorder eigenvalues, transform Generalized Schur vectors, and
674* compute reciprocal condition numbers
675*
676 CALL dtgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
677 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
678 $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
679 $ iwork, liwork, ierr )
680*
681 IF( ijob.GE.1 )
682 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
683 IF( ierr.EQ.-22 ) THEN
684*
685* not enough real workspace
686*
687 info = -22
688 ELSE
689 IF( ijob.EQ.1 .OR. ijob.EQ.4 ) THEN
690 rconde( 1 ) = pl
691 rconde( 2 ) = pr
692 END IF
693 IF( ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
694 rcondv( 1 ) = dif( 1 )
695 rcondv( 2 ) = dif( 2 )
696 END IF
697 IF( ierr.EQ.1 )
698 $ info = n + 3
699 END IF
700*
701 END IF
702*
703* Apply permutation to VSL and VSR
704* (Workspace: none needed)
705*
706 IF( ilvsl )
707 $ CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
708 $ work( iright ), n, vsl, ldvsl, ierr )
709*
710 IF( ilvsr )
711 $ CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
712 $ work( iright ), n, vsr, ldvsr, ierr )
713*
714* Check if unscaling would cause over/underflow, if so, rescale
715* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
716* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
717*
718 IF( ilascl ) THEN
719 DO 20 i = 1, n
720 IF( alphai( i ).NE.zero ) THEN
721 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
722 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) ) THEN
723 work( 1 ) = abs( a( i, i ) / alphar( i ) )
724 beta( i ) = beta( i )*work( 1 )
725 alphar( i ) = alphar( i )*work( 1 )
726 alphai( i ) = alphai( i )*work( 1 )
727 ELSE IF( ( alphai( i ) / safmax ).GT.
728 $ ( anrmto / anrm ) .OR.
729 $ ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
730 $ THEN
731 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
732 beta( i ) = beta( i )*work( 1 )
733 alphar( i ) = alphar( i )*work( 1 )
734 alphai( i ) = alphai( i )*work( 1 )
735 END IF
736 END IF
737 20 CONTINUE
738 END IF
739*
740 IF( ilbscl ) THEN
741 DO 30 i = 1, n
742 IF( alphai( i ).NE.zero ) THEN
743 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
744 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) ) THEN
745 work( 1 ) = abs( b( i, i ) / beta( i ) )
746 beta( i ) = beta( i )*work( 1 )
747 alphar( i ) = alphar( i )*work( 1 )
748 alphai( i ) = alphai( i )*work( 1 )
749 END IF
750 END IF
751 30 CONTINUE
752 END IF
753*
754* Undo scaling
755*
756 IF( ilascl ) THEN
757 CALL dlascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
758 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
759 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
760 END IF
761*
762 IF( ilbscl ) THEN
763 CALL dlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
764 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
765 END IF
766*
767 IF( wantst ) THEN
768*
769* Check if reordering is correct
770*
771 lastsl = .true.
772 lst2sl = .true.
773 sdim = 0
774 ip = 0
775 DO 50 i = 1, n
776 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
777 IF( alphai( i ).EQ.zero ) THEN
778 IF( cursl )
779 $ sdim = sdim + 1
780 ip = 0
781 IF( cursl .AND. .NOT.lastsl )
782 $ info = n + 2
783 ELSE
784 IF( ip.EQ.1 ) THEN
785*
786* Last eigenvalue of conjugate pair
787*
788 cursl = cursl .OR. lastsl
789 lastsl = cursl
790 IF( cursl )
791 $ sdim = sdim + 2
792 ip = -1
793 IF( cursl .AND. .NOT.lst2sl )
794 $ info = n + 2
795 ELSE
796*
797* First eigenvalue of conjugate pair
798*
799 ip = 1
800 END IF
801 END IF
802 lst2sl = lastsl
803 lastsl = cursl
804 50 CONTINUE
805*
806 END IF
807*
808 60 CONTINUE
809*
810 work( 1 ) = maxwrk
811 iwork( 1 ) = liwmin
812*
813 RETURN
814*
815* End of DGGESX
816*

◆ dggev()

subroutine dggev ( character jobvl,
character jobvr,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvl, * ) vl,
integer ldvl,
double precision, dimension( ldvr, * ) vr,
integer ldvr,
double precision, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> DGGEV computes for a pair of N-by-N real 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 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]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 DOUBLE PRECISION 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 DOUBLE PRECISION 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]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(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, ALPHAR and ALPHAI 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 DOUBLE PRECISION 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 the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
!>          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 DOUBLE PRECISION 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 the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
!>          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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,8*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]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 ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
!>                =N+2: error return from DTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file dggev.f.

226*
227* -- LAPACK driver routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232 CHARACTER JOBVL, JOBVR
233 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
234* ..
235* .. Array Arguments ..
236 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
237 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
238 $ VR( LDVR, * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 DOUBLE PRECISION ZERO, ONE
245 parameter( zero = 0.0d+0, one = 1.0d+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
249 CHARACTER CHTEMP
250 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
251 $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
252 $ MINWRK
253 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
254 $ SMLNUM, TEMP
255* ..
256* .. Local Arrays ..
257 LOGICAL LDUMMA( 1 )
258* ..
259* .. External Subroutines ..
260 EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
261 $ dlacpy,dlascl, dlaset, dorgqr, dormqr, dtgevc,
262 $ xerbla
263* ..
264* .. External Functions ..
265 LOGICAL LSAME
266 INTEGER ILAENV
267 DOUBLE PRECISION DLAMCH, DLANGE
268 EXTERNAL lsame, ilaenv, dlamch, dlange
269* ..
270* .. Intrinsic Functions ..
271 INTRINSIC abs, max, sqrt
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 = -12
316 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
317 info = -14
318 END IF
319*
320* Compute workspace
321* (Note: Comments in the code beginning "Workspace:" describe the
322* minimal amount of workspace needed at that point in the code,
323* as well as the preferred amount for good performance.
324* NB refers to the optimal block size for the immediately
325* following subroutine, as returned by ILAENV. The workspace is
326* computed assuming ILO = 1 and IHI = N, the worst case.)
327*
328 IF( info.EQ.0 ) THEN
329 minwrk = max( 1, 8*n )
330 maxwrk = max( 1, n*( 7 +
331 $ ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 ) ) )
332 maxwrk = max( maxwrk, n*( 7 +
333 $ ilaenv( 1, 'DORMQR', ' ', n, 1, n, 0 ) ) )
334 IF( ilvl ) THEN
335 maxwrk = max( maxwrk, n*( 7 +
336 $ ilaenv( 1, 'DORGQR', ' ', n, 1, n, -1 ) ) )
337 END IF
338 work( 1 ) = maxwrk
339*
340 IF( lwork.LT.minwrk .AND. .NOT.lquery )
341 $ info = -16
342 END IF
343*
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'DGGEV ', -info )
346 RETURN
347 ELSE IF( lquery ) THEN
348 RETURN
349 END IF
350*
351* Quick return if possible
352*
353 IF( n.EQ.0 )
354 $ RETURN
355*
356* Get machine constants
357*
358 eps = dlamch( 'P' )
359 smlnum = dlamch( 'S' )
360 bignum = one / smlnum
361 CALL dlabad( smlnum, bignum )
362 smlnum = sqrt( smlnum ) / eps
363 bignum = one / smlnum
364*
365* Scale A if max element outside range [SMLNUM,BIGNUM]
366*
367 anrm = dlange( 'M', n, n, a, lda, work )
368 ilascl = .false.
369 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
370 anrmto = smlnum
371 ilascl = .true.
372 ELSE IF( anrm.GT.bignum ) THEN
373 anrmto = bignum
374 ilascl = .true.
375 END IF
376 IF( ilascl )
377 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
378*
379* Scale B if max element outside range [SMLNUM,BIGNUM]
380*
381 bnrm = dlange( 'M', n, n, b, ldb, work )
382 ilbscl = .false.
383 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
384 bnrmto = smlnum
385 ilbscl = .true.
386 ELSE IF( bnrm.GT.bignum ) THEN
387 bnrmto = bignum
388 ilbscl = .true.
389 END IF
390 IF( ilbscl )
391 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
392*
393* Permute the matrices A, B to isolate eigenvalues if possible
394* (Workspace: need 6*N)
395*
396 ileft = 1
397 iright = n + 1
398 iwrk = iright + n
399 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
400 $ work( iright ), work( iwrk ), ierr )
401*
402* Reduce B to triangular form (QR decomposition of B)
403* (Workspace: need N, prefer N*NB)
404*
405 irows = ihi + 1 - ilo
406 IF( ilv ) THEN
407 icols = n + 1 - ilo
408 ELSE
409 icols = irows
410 END IF
411 itau = iwrk
412 iwrk = itau + irows
413 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
414 $ work( iwrk ), lwork+1-iwrk, ierr )
415*
416* Apply the orthogonal transformation to matrix A
417* (Workspace: need N, prefer N*NB)
418*
419 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
420 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
421 $ lwork+1-iwrk, ierr )
422*
423* Initialize VL
424* (Workspace: need N, prefer N*NB)
425*
426 IF( ilvl ) THEN
427 CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
428 IF( irows.GT.1 ) THEN
429 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
430 $ vl( ilo+1, ilo ), ldvl )
431 END IF
432 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
433 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
434 END IF
435*
436* Initialize VR
437*
438 IF( ilvr )
439 $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
440*
441* Reduce to generalized Hessenberg form
442* (Workspace: none needed)
443*
444 IF( ilv ) THEN
445*
446* Eigenvectors requested -- work on whole matrix.
447*
448 CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
449 $ ldvl, vr, ldvr, ierr )
450 ELSE
451 CALL dgghrd( '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 forms and Schur vectors)
457* (Workspace: need N)
458*
459 iwrk = itau
460 IF( ilv ) THEN
461 chtemp = 'S'
462 ELSE
463 chtemp = 'E'
464 END IF
465 CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
466 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
467 $ work( iwrk ), lwork+1-iwrk, ierr )
468 IF( ierr.NE.0 ) THEN
469 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
470 info = ierr
471 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
472 info = ierr - n
473 ELSE
474 info = n + 1
475 END IF
476 GO TO 110
477 END IF
478*
479* Compute Eigenvectors
480* (Workspace: need 6*N)
481*
482 IF( ilv ) THEN
483 IF( ilvl ) THEN
484 IF( ilvr ) THEN
485 chtemp = 'B'
486 ELSE
487 chtemp = 'L'
488 END IF
489 ELSE
490 chtemp = 'R'
491 END IF
492 CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
493 $ vr, ldvr, n, in, work( iwrk ), ierr )
494 IF( ierr.NE.0 ) THEN
495 info = n + 2
496 GO TO 110
497 END IF
498*
499* Undo balancing on VL and VR and normalization
500* (Workspace: none needed)
501*
502 IF( ilvl ) THEN
503 CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
504 $ work( iright ), n, vl, ldvl, ierr )
505 DO 50 jc = 1, n
506 IF( alphai( jc ).LT.zero )
507 $ GO TO 50
508 temp = zero
509 IF( alphai( jc ).EQ.zero ) THEN
510 DO 10 jr = 1, n
511 temp = max( temp, abs( vl( jr, jc ) ) )
512 10 CONTINUE
513 ELSE
514 DO 20 jr = 1, n
515 temp = max( temp, abs( vl( jr, jc ) )+
516 $ abs( vl( jr, jc+1 ) ) )
517 20 CONTINUE
518 END IF
519 IF( temp.LT.smlnum )
520 $ GO TO 50
521 temp = one / temp
522 IF( alphai( jc ).EQ.zero ) THEN
523 DO 30 jr = 1, n
524 vl( jr, jc ) = vl( jr, jc )*temp
525 30 CONTINUE
526 ELSE
527 DO 40 jr = 1, n
528 vl( jr, jc ) = vl( jr, jc )*temp
529 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
530 40 CONTINUE
531 END IF
532 50 CONTINUE
533 END IF
534 IF( ilvr ) THEN
535 CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
536 $ work( iright ), n, vr, ldvr, ierr )
537 DO 100 jc = 1, n
538 IF( alphai( jc ).LT.zero )
539 $ GO TO 100
540 temp = zero
541 IF( alphai( jc ).EQ.zero ) THEN
542 DO 60 jr = 1, n
543 temp = max( temp, abs( vr( jr, jc ) ) )
544 60 CONTINUE
545 ELSE
546 DO 70 jr = 1, n
547 temp = max( temp, abs( vr( jr, jc ) )+
548 $ abs( vr( jr, jc+1 ) ) )
549 70 CONTINUE
550 END IF
551 IF( temp.LT.smlnum )
552 $ GO TO 100
553 temp = one / temp
554 IF( alphai( jc ).EQ.zero ) THEN
555 DO 80 jr = 1, n
556 vr( jr, jc ) = vr( jr, jc )*temp
557 80 CONTINUE
558 ELSE
559 DO 90 jr = 1, n
560 vr( jr, jc ) = vr( jr, jc )*temp
561 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
562 90 CONTINUE
563 END IF
564 100 CONTINUE
565 END IF
566*
567* End of eigenvector calculation
568*
569 END IF
570*
571* Undo scaling if necessary
572*
573 110 CONTINUE
574*
575 IF( ilascl ) THEN
576 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
577 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
578 END IF
579*
580 IF( ilbscl ) THEN
581 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
582 END IF
583*
584 work( 1 ) = maxwrk
585 RETURN
586*
587* End of DGGEV
588*

◆ dggev3()

subroutine dggev3 ( character jobvl,
character jobvr,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvl, * ) vl,
integer ldvl,
double precision, dimension( ldvr, * ) vr,
integer ldvr,
double precision, dimension( * ) work,
integer lwork,
integer info )

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

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

Purpose:
!>
!> DGGEV3 computes for a pair of N-by-N real 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 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]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 DOUBLE PRECISION 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 DOUBLE PRECISION 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]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(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, ALPHAR and ALPHAI 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 DOUBLE PRECISION 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 the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
!>          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 DOUBLE PRECISION 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 the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
!>          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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>
!>          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]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 ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DLAQZ0.
!>                =N+2: error return from DTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file dggev3.f.

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

◆ dggevx()

subroutine dggevx ( character balanc,
character jobvl,
character jobvr,
character sense,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) alphar,
double precision, dimension( * ) alphai,
double precision, dimension( * ) beta,
double precision, dimension( ldvl, * ) vl,
integer ldvl,
double precision, dimension( ldvr, * ) vr,
integer ldvr,
integer ilo,
integer ihi,
double precision, dimension( * ) lscale,
double precision, dimension( * ) rscale,
double precision abnrm,
double precision bbnrm,
double precision, dimension( * ) rconde,
double precision, dimension( * ) rcondv,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
logical, dimension( * ) bwork,
integer info )

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

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

Purpose:
!>
!> DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
!> the generalized eigenvalues, and optionally, the left and/or right
!> generalized eigenvectors.
!>
!> Optionally also, it 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 DOUBLE PRECISION 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 real 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 DOUBLE PRECISION 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 real 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]ALPHAR
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!>
[out]BETA
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(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, ALPHAR and ALPHAI 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 DOUBLE PRECISION 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 the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
!>          Each eigenvector will be scaled so the largest component 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 DOUBLE PRECISION 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 the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
!>          Each eigenvector will be scaled so the largest component have
!>          abs(real part) + abs(imag. part) = 1.
!>          Not referenced if JOBVR = 'N'.
!>
[in]LDVR
!>          LDVR is INTEGER
!>          The leading dimension of the matrix VR. LDVR >= 1, and
!>          if JOBVR = 'V', LDVR >= N.
!>
[out]ILO
!>          ILO is INTEGER
!>
[out]IHI
!>          IHI is INTEGER
!>          ILO and IHI are integer values such that on exit
!>          A(i,j) = 0 and B(i,j) = 0 if i > j and
!>          j = 1,...,ILO-1 or i = IHI+1,...,N.
!>          If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
!>
[out]LSCALE
!>          LSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the left side of A and B.  If PL(j) is the index of the
!>          row interchanged with row j, and DL(j) is the scaling
!>          factor applied to row j, then
!>            LSCALE(j) = PL(j)  for j = 1,...,ILO-1
!>                      = DL(j)  for j = ILO,...,IHI
!>                      = PL(j)  for j = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]RSCALE
!>          RSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the right side of A and B.  If PR(j) is the index of the
!>          column interchanged with column j, and DR(j) is the scaling
!>          factor applied to column j, then
!>            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
!>                      = DR(j)  for j = ILO,...,IHI
!>                      = PR(j)  for j = IHI+1,...,N
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]ABNRM
!>          ABNRM is DOUBLE PRECISION
!>          The one-norm of the balanced matrix A.
!>
[out]BBNRM
!>          BBNRM is DOUBLE PRECISION
!>          The one-norm of the balanced matrix B.
!>
[out]RCONDE
!>          RCONDE is DOUBLE PRECISION array, dimension (N)
!>          If SENSE = 'E' or 'B', the reciprocal condition numbers of
!>          the eigenvalues, stored in consecutive elements of the array.
!>          For a complex conjugate pair of eigenvalues two consecutive
!>          elements of RCONDE are set to the same value. Thus RCONDE(j),
!>          RCONDV(j), and the j-th columns of VL and VR all correspond
!>          to the j-th eigenpair.
!>          If SENSE = 'N or 'V', RCONDE is not referenced.
!>
[out]RCONDV
!>          RCONDV is DOUBLE PRECISION array, dimension (N)
!>          If SENSE = 'V' or 'B', the estimated reciprocal condition
!>          numbers of the eigenvectors, stored in consecutive elements
!>          of the array. For a complex eigenvector two consecutive
!>          elements of RCONDV are set to the same value. 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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= max(1,2*N).
!>          If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
!>          LWORK >= max(1,6*N).
!>          If SENSE = 'E' or 'B', LWORK >= max(1,10*N).
!>          If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
!>
!>          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]IWORK
!>          IWORK is INTEGER array, dimension (N+6)
!>          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 ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
!>                =N+2: error return from DTGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 387 of file dggevx.f.

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