LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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,
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 = int( 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 smlnum = sqrt( smlnum ) / eps
342 bignum = one / smlnum
343*
344* Scale A if max element outside range [SMLNUM,BIGNUM]
345*
346 anrm = dlange( 'M', n, n, a, lda, dum )
347 scalea = .false.
348 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
349 scalea = .true.
350 cscale = smlnum
351 ELSE IF( anrm.GT.bignum ) THEN
352 scalea = .true.
353 cscale = bignum
354 END IF
355 IF( scalea )
356 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
357*
358* Permute the matrix to make it more nearly triangular
359* (Workspace: need N)
360*
361 ibal = 1
362 CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
363*
364* Reduce to upper Hessenberg form
365* (Workspace: need 3*N, prefer 2*N+N*NB)
366*
367 itau = n + ibal
368 iwrk = n + itau
369 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
370 $ lwork-iwrk+1, ierr )
371*
372 IF( wantvs ) THEN
373*
374* Copy Householder vectors to VS
375*
376 CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
377*
378* Generate orthogonal matrix in VS
379* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
380*
381 CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
382 $ lwork-iwrk+1, ierr )
383 END IF
384*
385 sdim = 0
386*
387* Perform QR iteration, accumulating Schur vectors in VS if desired
388* (Workspace: need N+1, prefer N+HSWORK (see comments) )
389*
390 iwrk = itau
391 CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
392 $ work( iwrk ), lwork-iwrk+1, ieval )
393 IF( ieval.GT.0 )
394 $ info = ieval
395*
396* Sort eigenvalues if desired
397*
398 IF( wantst .AND. info.EQ.0 ) THEN
399 IF( scalea ) THEN
400 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
401 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
402 END IF
403 DO 10 i = 1, n
404 bwork( i ) = SELECT( wr( i ), wi( i ) )
405 10 CONTINUE
406*
407* Reorder eigenvalues and transform Schur vectors
408* (Workspace: none needed)
409*
410 CALL dtrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
411 $ sdim, s, sep, work( iwrk ), lwork-iwrk+1, idum, 1,
412 $ icond )
413 IF( icond.GT.0 )
414 $ info = n + icond
415 END IF
416*
417 IF( wantvs ) THEN
418*
419* Undo balancing
420* (Workspace: need N)
421*
422 CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
423 $ ierr )
424 END IF
425*
426 IF( scalea ) THEN
427*
428* Undo scaling for the Schur form of A
429*
430 CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
431 CALL dcopy( n, a, lda+1, wr, 1 )
432 IF( cscale.EQ.smlnum ) THEN
433*
434* If scaling back towards underflow, adjust WI if an
435* offdiagonal element of a 2-by-2 block in the Schur form
436* underflows.
437*
438 IF( ieval.GT.0 ) THEN
439 i1 = ieval + 1
440 i2 = ihi - 1
441 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi,
442 $ max( ilo-1, 1 ), ierr )
443 ELSE IF( wantst ) THEN
444 i1 = 1
445 i2 = n - 1
446 ELSE
447 i1 = ilo
448 i2 = ihi - 1
449 END IF
450 inxt = i1 - 1
451 DO 20 i = i1, i2
452 IF( i.LT.inxt )
453 $ GO TO 20
454 IF( wi( i ).EQ.zero ) THEN
455 inxt = i + 1
456 ELSE
457 IF( a( i+1, i ).EQ.zero ) THEN
458 wi( i ) = zero
459 wi( i+1 ) = zero
460 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
461 $ zero ) THEN
462 wi( i ) = zero
463 wi( i+1 ) = zero
464 IF( i.GT.1 )
465 $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
466 IF( n.GT.i+1 )
467 $ CALL dswap( n-i-1, a( i, i+2 ), lda,
468 $ a( i+1, i+2 ), lda )
469 IF( wantvs ) THEN
470 CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
471 END IF
472 a( i, i+1 ) = a( i+1, i )
473 a( i+1, i ) = zero
474 END IF
475 inxt = i + 2
476 END IF
477 20 CONTINUE
478 END IF
479*
480* Undo scaling for the imaginary part of the eigenvalues
481*
482 CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
483 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
484 END IF
485*
486 IF( wantst .AND. info.EQ.0 ) THEN
487*
488* Check if reordering successful
489*
490 lastsl = .true.
491 lst2sl = .true.
492 sdim = 0
493 ip = 0
494 DO 30 i = 1, n
495 cursl = SELECT( wr( i ), wi( i ) )
496 IF( wi( i ).EQ.zero ) THEN
497 IF( cursl )
498 $ sdim = sdim + 1
499 ip = 0
500 IF( cursl .AND. .NOT.lastsl )
501 $ info = n + 2
502 ELSE
503 IF( ip.EQ.1 ) THEN
504*
505* Last eigenvalue of conjugate pair
506*
507 cursl = cursl .OR. lastsl
508 lastsl = cursl
509 IF( cursl )
510 $ sdim = sdim + 2
511 ip = -1
512 IF( cursl .AND. .NOT.lst2sl )
513 $ info = n + 2
514 ELSE
515*
516* First eigenvalue of conjugate pair
517*
518 ip = 1
519 END IF
520 END IF
521 lst2sl = lastsl
522 lastsl = cursl
523 30 CONTINUE
524 END IF
525*
526 work( 1 ) = maxwrk
527 RETURN
528*
529* End of DGEES
530*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
DGEBAK
Definition dgebak.f:130
subroutine dgebal(job, n, a, lda, ilo, ihi, scale, info)
DGEBAL
Definition dgebal.f:163
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:316
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
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
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: