Annotation of rpl/lapack/lapack/dgees.f, revision 1.8
1.8 ! bertrand 1: *> \brief <b> DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGEES + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgees.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgees.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgees.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
! 22: * VS, LDVS, WORK, LWORK, BWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER JOBVS, SORT
! 26: * INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * LOGICAL BWORK( * )
! 30: * DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
! 31: * $ WR( * )
! 32: * ..
! 33: * .. Function Arguments ..
! 34: * LOGICAL SELECT
! 35: * EXTERNAL SELECT
! 36: * ..
! 37: *
! 38: *
! 39: *> \par Purpose:
! 40: * =============
! 41: *>
! 42: *> \verbatim
! 43: *>
! 44: *> DGEES computes for an N-by-N real nonsymmetric matrix A, the
! 45: *> eigenvalues, the real Schur form T, and, optionally, the matrix of
! 46: *> Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
! 47: *>
! 48: *> Optionally, it also orders the eigenvalues on the diagonal of the
! 49: *> real Schur form so that selected eigenvalues are at the top left.
! 50: *> The leading columns of Z then form an orthonormal basis for the
! 51: *> invariant subspace corresponding to the selected eigenvalues.
! 52: *>
! 53: *> A matrix is in real Schur form if it is upper quasi-triangular with
! 54: *> 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
! 55: *> form
! 56: *> [ a b ]
! 57: *> [ c a ]
! 58: *>
! 59: *> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
! 60: *> \endverbatim
! 61: *
! 62: * Arguments:
! 63: * ==========
! 64: *
! 65: *> \param[in] JOBVS
! 66: *> \verbatim
! 67: *> JOBVS is CHARACTER*1
! 68: *> = 'N': Schur vectors are not computed;
! 69: *> = 'V': Schur vectors are computed.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] SORT
! 73: *> \verbatim
! 74: *> SORT is CHARACTER*1
! 75: *> Specifies whether or not to order the eigenvalues on the
! 76: *> diagonal of the Schur form.
! 77: *> = 'N': Eigenvalues are not ordered;
! 78: *> = 'S': Eigenvalues are ordered (see SELECT).
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] SELECT
! 82: *> \verbatim
! 83: *> SELECT is procedure) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
! 84: *> SELECT must be declared EXTERNAL in the calling subroutine.
! 85: *> If SORT = 'S', SELECT is used to select eigenvalues to sort
! 86: *> to the top left of the Schur form.
! 87: *> If SORT = 'N', SELECT is not referenced.
! 88: *> An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
! 89: *> SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
! 90: *> conjugate pair of eigenvalues is selected, then both complex
! 91: *> eigenvalues are selected.
! 92: *> Note that a selected complex eigenvalue may no longer
! 93: *> satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
! 94: *> ordering may change the value of complex eigenvalues
! 95: *> (especially if the eigenvalue is ill-conditioned); in this
! 96: *> case INFO is set to N+2 (see INFO below).
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] N
! 100: *> \verbatim
! 101: *> N is INTEGER
! 102: *> The order of the matrix A. N >= 0.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in,out] A
! 106: *> \verbatim
! 107: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 108: *> On entry, the N-by-N matrix A.
! 109: *> On exit, A has been overwritten by its real Schur form T.
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[in] LDA
! 113: *> \verbatim
! 114: *> LDA is INTEGER
! 115: *> The leading dimension of the array A. LDA >= max(1,N).
! 116: *> \endverbatim
! 117: *>
! 118: *> \param[out] SDIM
! 119: *> \verbatim
! 120: *> SDIM is INTEGER
! 121: *> If SORT = 'N', SDIM = 0.
! 122: *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
! 123: *> for which SELECT is true. (Complex conjugate
! 124: *> pairs for which SELECT is true for either
! 125: *> eigenvalue count as 2.)
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[out] WR
! 129: *> \verbatim
! 130: *> WR is DOUBLE PRECISION array, dimension (N)
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[out] WI
! 134: *> \verbatim
! 135: *> WI is DOUBLE PRECISION array, dimension (N)
! 136: *> WR and WI contain the real and imaginary parts,
! 137: *> respectively, of the computed eigenvalues in the same order
! 138: *> that they appear on the diagonal of the output Schur form T.
! 139: *> Complex conjugate pairs of eigenvalues will appear
! 140: *> consecutively with the eigenvalue having the positive
! 141: *> imaginary part first.
! 142: *> \endverbatim
! 143: *>
! 144: *> \param[out] VS
! 145: *> \verbatim
! 146: *> VS is DOUBLE PRECISION array, dimension (LDVS,N)
! 147: *> If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
! 148: *> vectors.
! 149: *> If JOBVS = 'N', VS is not referenced.
! 150: *> \endverbatim
! 151: *>
! 152: *> \param[in] LDVS
! 153: *> \verbatim
! 154: *> LDVS is INTEGER
! 155: *> The leading dimension of the array VS. LDVS >= 1; if
! 156: *> JOBVS = 'V', LDVS >= N.
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[out] WORK
! 160: *> \verbatim
! 161: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 162: *> On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
! 163: *> \endverbatim
! 164: *>
! 165: *> \param[in] LWORK
! 166: *> \verbatim
! 167: *> LWORK is INTEGER
! 168: *> The dimension of the array WORK. LWORK >= max(1,3*N).
! 169: *> For good performance, LWORK must generally be larger.
! 170: *>
! 171: *> If LWORK = -1, then a workspace query is assumed; the routine
! 172: *> only calculates the optimal size of the WORK array, returns
! 173: *> this value as the first entry of the WORK array, and no error
! 174: *> message related to LWORK is issued by XERBLA.
! 175: *> \endverbatim
! 176: *>
! 177: *> \param[out] BWORK
! 178: *> \verbatim
! 179: *> BWORK is LOGICAL array, dimension (N)
! 180: *> Not referenced if SORT = 'N'.
! 181: *> \endverbatim
! 182: *>
! 183: *> \param[out] INFO
! 184: *> \verbatim
! 185: *> INFO is INTEGER
! 186: *> = 0: successful exit
! 187: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 188: *> > 0: if INFO = i, and i is
! 189: *> <= N: the QR algorithm failed to compute all the
! 190: *> eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
! 191: *> contain those eigenvalues which have converged; if
! 192: *> JOBVS = 'V', VS contains the matrix which reduces A
! 193: *> to its partially converged Schur form.
! 194: *> = N+1: the eigenvalues could not be reordered because some
! 195: *> eigenvalues were too close to separate (the problem
! 196: *> is very ill-conditioned);
! 197: *> = N+2: after reordering, roundoff changed values of some
! 198: *> complex eigenvalues so that leading eigenvalues in
! 199: *> the Schur form no longer satisfy SELECT=.TRUE. This
! 200: *> could also be caused by underflow due to scaling.
! 201: *> \endverbatim
! 202: *
! 203: * Authors:
! 204: * ========
! 205: *
! 206: *> \author Univ. of Tennessee
! 207: *> \author Univ. of California Berkeley
! 208: *> \author Univ. of Colorado Denver
! 209: *> \author NAG Ltd.
! 210: *
! 211: *> \date November 2011
! 212: *
! 213: *> \ingroup doubleGEeigen
! 214: *
! 215: * =====================================================================
1.1 bertrand 216: SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
217: $ VS, LDVS, WORK, LWORK, BWORK, INFO )
218: *
1.8 ! bertrand 219: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 220: * -- LAPACK is a software package provided by Univ. of Tennessee, --
221: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 222: * November 2011
1.1 bertrand 223: *
224: * .. Scalar Arguments ..
225: CHARACTER JOBVS, SORT
226: INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
227: * ..
228: * .. Array Arguments ..
229: LOGICAL BWORK( * )
230: DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
231: $ WR( * )
232: * ..
233: * .. Function Arguments ..
234: LOGICAL SELECT
235: EXTERNAL SELECT
236: * ..
237: *
238: * =====================================================================
239: *
240: * .. Parameters ..
241: DOUBLE PRECISION ZERO, ONE
242: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
243: * ..
244: * .. Local Scalars ..
245: LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
246: $ WANTVS
247: INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
248: $ IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
249: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
250: * ..
251: * .. Local Arrays ..
252: INTEGER IDUM( 1 )
253: DOUBLE PRECISION DUM( 1 )
254: * ..
255: * .. External Subroutines ..
256: EXTERNAL DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY,
257: $ DLABAD, DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
258: * ..
259: * .. External Functions ..
260: LOGICAL LSAME
261: INTEGER ILAENV
262: DOUBLE PRECISION DLAMCH, DLANGE
263: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
264: * ..
265: * .. Intrinsic Functions ..
266: INTRINSIC MAX, SQRT
267: * ..
268: * .. Executable Statements ..
269: *
270: * Test the input arguments
271: *
272: INFO = 0
273: LQUERY = ( LWORK.EQ.-1 )
274: WANTVS = LSAME( JOBVS, 'V' )
275: WANTST = LSAME( SORT, 'S' )
276: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
277: INFO = -1
278: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
279: INFO = -2
280: ELSE IF( N.LT.0 ) THEN
281: INFO = -4
282: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
283: INFO = -6
284: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
285: INFO = -11
286: END IF
287: *
288: * Compute workspace
289: * (Note: Comments in the code beginning "Workspace:" describe the
290: * minimal amount of workspace needed at that point in the code,
291: * as well as the preferred amount for good performance.
292: * NB refers to the optimal block size for the immediately
293: * following subroutine, as returned by ILAENV.
294: * HSWORK refers to the workspace preferred by DHSEQR, as
295: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
296: * the worst case.)
297: *
298: IF( INFO.EQ.0 ) THEN
299: IF( N.EQ.0 ) THEN
300: MINWRK = 1
301: MAXWRK = 1
302: ELSE
303: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
304: MINWRK = 3*N
305: *
306: CALL DHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
307: $ WORK, -1, IEVAL )
308: HSWORK = WORK( 1 )
309: *
310: IF( .NOT.WANTVS ) THEN
311: MAXWRK = MAX( MAXWRK, N + HSWORK )
312: ELSE
313: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
314: $ 'DORGHR', ' ', N, 1, N, -1 ) )
315: MAXWRK = MAX( MAXWRK, N + HSWORK )
316: END IF
317: END IF
318: WORK( 1 ) = MAXWRK
319: *
320: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
321: INFO = -13
322: END IF
323: END IF
324: *
325: IF( INFO.NE.0 ) THEN
326: CALL XERBLA( 'DGEES ', -INFO )
327: RETURN
328: ELSE IF( LQUERY ) THEN
329: RETURN
330: END IF
331: *
332: * Quick return if possible
333: *
334: IF( N.EQ.0 ) THEN
335: SDIM = 0
336: RETURN
337: END IF
338: *
339: * Get machine constants
340: *
341: EPS = DLAMCH( 'P' )
342: SMLNUM = DLAMCH( 'S' )
343: BIGNUM = ONE / SMLNUM
344: CALL DLABAD( SMLNUM, BIGNUM )
345: SMLNUM = SQRT( SMLNUM ) / 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, DUM )
351: SCALEA = .FALSE.
352: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
353: SCALEA = .TRUE.
354: CSCALE = SMLNUM
355: ELSE IF( ANRM.GT.BIGNUM ) THEN
356: SCALEA = .TRUE.
357: CSCALE = BIGNUM
358: END IF
359: IF( SCALEA )
360: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
361: *
362: * Permute the matrix to make it more nearly triangular
363: * (Workspace: need N)
364: *
365: IBAL = 1
366: CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
367: *
368: * Reduce to upper Hessenberg form
369: * (Workspace: need 3*N, prefer 2*N+N*NB)
370: *
371: ITAU = N + IBAL
372: IWRK = N + ITAU
373: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
374: $ LWORK-IWRK+1, IERR )
375: *
376: IF( WANTVS ) THEN
377: *
378: * Copy Householder vectors to VS
379: *
380: CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
381: *
382: * Generate orthogonal matrix in VS
383: * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
384: *
385: CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
386: $ LWORK-IWRK+1, IERR )
387: END IF
388: *
389: SDIM = 0
390: *
391: * Perform QR iteration, accumulating Schur vectors in VS if desired
392: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
393: *
394: IWRK = ITAU
395: CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
396: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
397: IF( IEVAL.GT.0 )
398: $ INFO = IEVAL
399: *
400: * Sort eigenvalues if desired
401: *
402: IF( WANTST .AND. INFO.EQ.0 ) THEN
403: IF( SCALEA ) THEN
404: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
405: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
406: END IF
407: DO 10 I = 1, N
408: BWORK( I ) = SELECT( WR( I ), WI( I ) )
409: 10 CONTINUE
410: *
411: * Reorder eigenvalues and transform Schur vectors
412: * (Workspace: none needed)
413: *
414: CALL DTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
415: $ SDIM, S, SEP, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
416: $ ICOND )
417: IF( ICOND.GT.0 )
418: $ INFO = N + ICOND
419: END IF
420: *
421: IF( WANTVS ) THEN
422: *
423: * Undo balancing
424: * (Workspace: need N)
425: *
426: CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
427: $ IERR )
428: END IF
429: *
430: IF( SCALEA ) THEN
431: *
432: * Undo scaling for the Schur form of A
433: *
434: CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
435: CALL DCOPY( N, A, LDA+1, WR, 1 )
436: IF( CSCALE.EQ.SMLNUM ) THEN
437: *
438: * If scaling back towards underflow, adjust WI if an
439: * offdiagonal element of a 2-by-2 block in the Schur form
440: * underflows.
441: *
442: IF( IEVAL.GT.0 ) THEN
443: I1 = IEVAL + 1
444: I2 = IHI - 1
445: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI,
446: $ MAX( ILO-1, 1 ), IERR )
447: ELSE IF( WANTST ) THEN
448: I1 = 1
449: I2 = N - 1
450: ELSE
451: I1 = ILO
452: I2 = IHI - 1
453: END IF
454: INXT = I1 - 1
455: DO 20 I = I1, I2
456: IF( I.LT.INXT )
457: $ GO TO 20
458: IF( WI( I ).EQ.ZERO ) THEN
459: INXT = I + 1
460: ELSE
461: IF( A( I+1, I ).EQ.ZERO ) THEN
462: WI( I ) = ZERO
463: WI( I+1 ) = ZERO
464: ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
465: $ ZERO ) THEN
466: WI( I ) = ZERO
467: WI( I+1 ) = ZERO
468: IF( I.GT.1 )
469: $ CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
470: IF( N.GT.I+1 )
471: $ CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
472: $ A( I+1, I+2 ), LDA )
473: IF( WANTVS ) THEN
474: CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
475: END IF
476: A( I, I+1 ) = A( I+1, I )
477: A( I+1, I ) = ZERO
478: END IF
479: INXT = I + 2
480: END IF
481: 20 CONTINUE
482: END IF
483: *
484: * Undo scaling for the imaginary part of the eigenvalues
485: *
486: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
487: $ WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
488: END IF
489: *
490: IF( WANTST .AND. INFO.EQ.0 ) THEN
491: *
492: * Check if reordering successful
493: *
494: LASTSL = .TRUE.
495: LST2SL = .TRUE.
496: SDIM = 0
497: IP = 0
498: DO 30 I = 1, N
499: CURSL = SELECT( WR( I ), WI( I ) )
500: IF( WI( I ).EQ.ZERO ) THEN
501: IF( CURSL )
502: $ SDIM = SDIM + 1
503: IP = 0
504: IF( CURSL .AND. .NOT.LASTSL )
505: $ INFO = N + 2
506: ELSE
507: IF( IP.EQ.1 ) THEN
508: *
509: * Last eigenvalue of conjugate pair
510: *
511: CURSL = CURSL .OR. LASTSL
512: LASTSL = CURSL
513: IF( CURSL )
514: $ SDIM = SDIM + 2
515: IP = -1
516: IF( CURSL .AND. .NOT.LST2SL )
517: $ INFO = N + 2
518: ELSE
519: *
520: * First eigenvalue of conjugate pair
521: *
522: IP = 1
523: END IF
524: END IF
525: LST2SL = LASTSL
526: LASTSL = CURSL
527: 30 CONTINUE
528: END IF
529: *
530: WORK( 1 ) = MAXWRK
531: RETURN
532: *
533: * End of DGEES
534: *
535: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgees.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack