Annotation of rpl/lapack/lapack/zgees.f, revision 1.8
1.8 ! bertrand 1: *> \brief <b> ZGEES 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 ZGEES + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgees.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgees.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgees.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
! 22: * LDVS, WORK, LWORK, RWORK, 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 RWORK( * )
! 31: * COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
! 32: * ..
! 33: * .. Function Arguments ..
! 34: * LOGICAL SELECT
! 35: * EXTERNAL SELECT
! 36: * ..
! 37: *
! 38: *
! 39: *> \par Purpose:
! 40: * =============
! 41: *>
! 42: *> \verbatim
! 43: *>
! 44: *> ZGEES computes for an N-by-N complex nonsymmetric matrix A, the
! 45: *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
! 46: *> vectors Z. This gives the Schur factorization A = Z*T*(Z**H).
! 47: *>
! 48: *> Optionally, it also orders the eigenvalues on the diagonal of the
! 49: *> 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 complex matrix is in Schur form if it is upper triangular.
! 54: *> \endverbatim
! 55: *
! 56: * Arguments:
! 57: * ==========
! 58: *
! 59: *> \param[in] JOBVS
! 60: *> \verbatim
! 61: *> JOBVS is CHARACTER*1
! 62: *> = 'N': Schur vectors are not computed;
! 63: *> = 'V': Schur vectors are computed.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] SORT
! 67: *> \verbatim
! 68: *> SORT is CHARACTER*1
! 69: *> Specifies whether or not to order the eigenvalues on the
! 70: *> diagonal of the Schur form.
! 71: *> = 'N': Eigenvalues are not ordered:
! 72: *> = 'S': Eigenvalues are ordered (see SELECT).
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in] SELECT
! 76: *> \verbatim
! 77: *> SELECT is procedure) LOGICAL FUNCTION of one COMPLEX*16 argument
! 78: *> SELECT must be declared EXTERNAL in the calling subroutine.
! 79: *> If SORT = 'S', SELECT is used to select eigenvalues to order
! 80: *> to the top left of the Schur form.
! 81: *> IF SORT = 'N', SELECT is not referenced.
! 82: *> The eigenvalue W(j) is selected if SELECT(W(j)) is true.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in] N
! 86: *> \verbatim
! 87: *> N is INTEGER
! 88: *> The order of the matrix A. N >= 0.
! 89: *> \endverbatim
! 90: *>
! 91: *> \param[in,out] A
! 92: *> \verbatim
! 93: *> A is COMPLEX*16 array, dimension (LDA,N)
! 94: *> On entry, the N-by-N matrix A.
! 95: *> On exit, A has been overwritten by its Schur form T.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in] LDA
! 99: *> \verbatim
! 100: *> LDA is INTEGER
! 101: *> The leading dimension of the array A. LDA >= max(1,N).
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[out] SDIM
! 105: *> \verbatim
! 106: *> SDIM is INTEGER
! 107: *> If SORT = 'N', SDIM = 0.
! 108: *> If SORT = 'S', SDIM = number of eigenvalues for which
! 109: *> SELECT is true.
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[out] W
! 113: *> \verbatim
! 114: *> W is COMPLEX*16 array, dimension (N)
! 115: *> W contains the computed eigenvalues, in the same order that
! 116: *> they appear on the diagonal of the output Schur form T.
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[out] VS
! 120: *> \verbatim
! 121: *> VS is COMPLEX*16 array, dimension (LDVS,N)
! 122: *> If JOBVS = 'V', VS contains the unitary matrix Z of Schur
! 123: *> vectors.
! 124: *> If JOBVS = 'N', VS is not referenced.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in] LDVS
! 128: *> \verbatim
! 129: *> LDVS is INTEGER
! 130: *> The leading dimension of the array VS. LDVS >= 1; if
! 131: *> JOBVS = 'V', LDVS >= N.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[out] WORK
! 135: *> \verbatim
! 136: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 137: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[in] LWORK
! 141: *> \verbatim
! 142: *> LWORK is INTEGER
! 143: *> The dimension of the array WORK. LWORK >= max(1,2*N).
! 144: *> For good performance, LWORK must generally be larger.
! 145: *>
! 146: *> If LWORK = -1, then a workspace query is assumed; the routine
! 147: *> only calculates the optimal size of the WORK array, returns
! 148: *> this value as the first entry of the WORK array, and no error
! 149: *> message related to LWORK is issued by XERBLA.
! 150: *> \endverbatim
! 151: *>
! 152: *> \param[out] RWORK
! 153: *> \verbatim
! 154: *> RWORK is DOUBLE PRECISION array, dimension (N)
! 155: *> \endverbatim
! 156: *>
! 157: *> \param[out] BWORK
! 158: *> \verbatim
! 159: *> BWORK is LOGICAL array, dimension (N)
! 160: *> Not referenced if SORT = 'N'.
! 161: *> \endverbatim
! 162: *>
! 163: *> \param[out] INFO
! 164: *> \verbatim
! 165: *> INFO is INTEGER
! 166: *> = 0: successful exit
! 167: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 168: *> > 0: if INFO = i, and i is
! 169: *> <= N: the QR algorithm failed to compute all the
! 170: *> eigenvalues; elements 1:ILO-1 and i+1:N of W
! 171: *> contain those eigenvalues which have converged;
! 172: *> if JOBVS = 'V', VS contains the matrix which
! 173: *> reduces A to its partially converged Schur form.
! 174: *> = N+1: the eigenvalues could not be reordered because
! 175: *> some eigenvalues were too close to separate (the
! 176: *> problem is very ill-conditioned);
! 177: *> = N+2: after reordering, roundoff changed values of
! 178: *> some complex eigenvalues so that leading
! 179: *> eigenvalues in the Schur form no longer satisfy
! 180: *> SELECT = .TRUE.. This could also be caused by
! 181: *> underflow due to scaling.
! 182: *> \endverbatim
! 183: *
! 184: * Authors:
! 185: * ========
! 186: *
! 187: *> \author Univ. of Tennessee
! 188: *> \author Univ. of California Berkeley
! 189: *> \author Univ. of Colorado Denver
! 190: *> \author NAG Ltd.
! 191: *
! 192: *> \date November 2011
! 193: *
! 194: *> \ingroup complex16GEeigen
! 195: *
! 196: * =====================================================================
1.1 bertrand 197: SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
198: $ LDVS, WORK, LWORK, RWORK, BWORK, INFO )
199: *
1.8 ! bertrand 200: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 201: * -- LAPACK is a software package provided by Univ. of Tennessee, --
202: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 203: * November 2011
1.1 bertrand 204: *
205: * .. Scalar Arguments ..
206: CHARACTER JOBVS, SORT
207: INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
208: * ..
209: * .. Array Arguments ..
210: LOGICAL BWORK( * )
211: DOUBLE PRECISION RWORK( * )
212: COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
213: * ..
214: * .. Function Arguments ..
215: LOGICAL SELECT
216: EXTERNAL SELECT
217: * ..
218: *
219: * =====================================================================
220: *
221: * .. Parameters ..
222: DOUBLE PRECISION ZERO, ONE
223: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
224: * ..
225: * .. Local Scalars ..
226: LOGICAL LQUERY, SCALEA, WANTST, WANTVS
227: INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
228: $ ITAU, IWRK, MAXWRK, MINWRK
229: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
230: * ..
231: * .. Local Arrays ..
232: DOUBLE PRECISION DUM( 1 )
233: * ..
234: * .. External Subroutines ..
235: EXTERNAL DLABAD, XERBLA, ZCOPY, ZGEBAK, ZGEBAL, ZGEHRD,
236: $ ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
237: * ..
238: * .. External Functions ..
239: LOGICAL LSAME
240: INTEGER ILAENV
241: DOUBLE PRECISION DLAMCH, ZLANGE
242: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
243: * ..
244: * .. Intrinsic Functions ..
245: INTRINSIC MAX, SQRT
246: * ..
247: * .. Executable Statements ..
248: *
249: * Test the input arguments
250: *
251: INFO = 0
252: LQUERY = ( LWORK.EQ.-1 )
253: WANTVS = LSAME( JOBVS, 'V' )
254: WANTST = LSAME( SORT, 'S' )
255: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
256: INFO = -1
257: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
258: INFO = -2
259: ELSE IF( N.LT.0 ) THEN
260: INFO = -4
261: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
262: INFO = -6
263: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
264: INFO = -10
265: END IF
266: *
267: * Compute workspace
268: * (Note: Comments in the code beginning "Workspace:" describe the
269: * minimal amount of workspace needed at that point in the code,
270: * as well as the preferred amount for good performance.
271: * CWorkspace refers to complex workspace, and RWorkspace to real
272: * workspace. NB refers to the optimal block size for the
273: * immediately following subroutine, as returned by ILAENV.
274: * HSWORK refers to the workspace preferred by ZHSEQR, as
275: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
276: * the worst case.)
277: *
278: IF( INFO.EQ.0 ) THEN
279: IF( N.EQ.0 ) THEN
280: MINWRK = 1
281: MAXWRK = 1
282: ELSE
283: MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
284: MINWRK = 2*N
285: *
286: CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
287: $ WORK, -1, IEVAL )
288: HSWORK = WORK( 1 )
289: *
290: IF( .NOT.WANTVS ) THEN
291: MAXWRK = MAX( MAXWRK, HSWORK )
292: ELSE
293: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
294: $ ' ', N, 1, N, -1 ) )
295: MAXWRK = MAX( MAXWRK, HSWORK )
296: END IF
297: END IF
298: WORK( 1 ) = MAXWRK
299: *
300: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
301: INFO = -12
302: END IF
303: END IF
304: *
305: IF( INFO.NE.0 ) THEN
306: CALL XERBLA( 'ZGEES ', -INFO )
307: RETURN
308: ELSE IF( LQUERY ) THEN
309: RETURN
310: END IF
311: *
312: * Quick return if possible
313: *
314: IF( N.EQ.0 ) THEN
315: SDIM = 0
316: RETURN
317: END IF
318: *
319: * Get machine constants
320: *
321: EPS = DLAMCH( 'P' )
322: SMLNUM = DLAMCH( 'S' )
323: BIGNUM = ONE / SMLNUM
324: CALL DLABAD( SMLNUM, BIGNUM )
325: SMLNUM = SQRT( SMLNUM ) / EPS
326: BIGNUM = ONE / SMLNUM
327: *
328: * Scale A if max element outside range [SMLNUM,BIGNUM]
329: *
330: ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
331: SCALEA = .FALSE.
332: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
333: SCALEA = .TRUE.
334: CSCALE = SMLNUM
335: ELSE IF( ANRM.GT.BIGNUM ) THEN
336: SCALEA = .TRUE.
337: CSCALE = BIGNUM
338: END IF
339: IF( SCALEA )
340: $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
341: *
342: * Permute the matrix to make it more nearly triangular
343: * (CWorkspace: none)
344: * (RWorkspace: need N)
345: *
346: IBAL = 1
347: CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
348: *
349: * Reduce to upper Hessenberg form
350: * (CWorkspace: need 2*N, prefer N+N*NB)
351: * (RWorkspace: none)
352: *
353: ITAU = 1
354: IWRK = N + ITAU
355: CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
356: $ LWORK-IWRK+1, IERR )
357: *
358: IF( WANTVS ) THEN
359: *
360: * Copy Householder vectors to VS
361: *
362: CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
363: *
364: * Generate unitary matrix in VS
365: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
366: * (RWorkspace: none)
367: *
368: CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
369: $ LWORK-IWRK+1, IERR )
370: END IF
371: *
372: SDIM = 0
373: *
374: * Perform QR iteration, accumulating Schur vectors in VS if desired
375: * (CWorkspace: need 1, prefer HSWORK (see comments) )
376: * (RWorkspace: none)
377: *
378: IWRK = ITAU
379: CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
380: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
381: IF( IEVAL.GT.0 )
382: $ INFO = IEVAL
383: *
384: * Sort eigenvalues if desired
385: *
386: IF( WANTST .AND. INFO.EQ.0 ) THEN
387: IF( SCALEA )
388: $ CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
389: DO 10 I = 1, N
390: BWORK( I ) = SELECT( W( I ) )
391: 10 CONTINUE
392: *
393: * Reorder eigenvalues and transform Schur vectors
394: * (CWorkspace: none)
395: * (RWorkspace: none)
396: *
397: CALL ZTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
398: $ S, SEP, WORK( IWRK ), LWORK-IWRK+1, ICOND )
399: END IF
400: *
401: IF( WANTVS ) THEN
402: *
403: * Undo balancing
404: * (CWorkspace: none)
405: * (RWorkspace: need N)
406: *
407: CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
408: $ IERR )
409: END IF
410: *
411: IF( SCALEA ) THEN
412: *
413: * Undo scaling for the Schur form of A
414: *
415: CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
416: CALL ZCOPY( N, A, LDA+1, W, 1 )
417: END IF
418: *
419: WORK( 1 ) = MAXWRK
420: RETURN
421: *
422: * End of ZGEES
423: *
424: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgees.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack