Annotation of rpl/lapack/lapack/dsygvd.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DSYGST
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DSYGVD + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygvd.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygvd.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygvd.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
! 22: * LWORK, IWORK, LIWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER JOBZ, UPLO
! 26: * INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IWORK( * )
! 30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DSYGVD computes all the eigenvalues, and optionally, the eigenvectors
! 40: *> of a real generalized symmetric-definite eigenproblem, of the form
! 41: *> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
! 42: *> B are assumed to be symmetric and B is also positive definite.
! 43: *> If eigenvectors are desired, it uses a divide and conquer algorithm.
! 44: *>
! 45: *> The divide and conquer algorithm makes very mild assumptions about
! 46: *> floating point arithmetic. It will work on machines with a guard
! 47: *> digit in add/subtract, or on those binary machines without guard
! 48: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
! 49: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
! 50: *> without guard digits, but we know of none.
! 51: *> \endverbatim
! 52: *
! 53: * Arguments:
! 54: * ==========
! 55: *
! 56: *> \param[in] ITYPE
! 57: *> \verbatim
! 58: *> ITYPE is INTEGER
! 59: *> Specifies the problem type to be solved:
! 60: *> = 1: A*x = (lambda)*B*x
! 61: *> = 2: A*B*x = (lambda)*x
! 62: *> = 3: B*A*x = (lambda)*x
! 63: *> \endverbatim
! 64: *>
! 65: *> \param[in] JOBZ
! 66: *> \verbatim
! 67: *> JOBZ is CHARACTER*1
! 68: *> = 'N': Compute eigenvalues only;
! 69: *> = 'V': Compute eigenvalues and eigenvectors.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] UPLO
! 73: *> \verbatim
! 74: *> UPLO is CHARACTER*1
! 75: *> = 'U': Upper triangles of A and B are stored;
! 76: *> = 'L': Lower triangles of A and B are stored.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] N
! 80: *> \verbatim
! 81: *> N is INTEGER
! 82: *> The order of the matrices A and B. N >= 0.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in,out] A
! 86: *> \verbatim
! 87: *> A is DOUBLE PRECISION array, dimension (LDA, N)
! 88: *> On entry, the symmetric matrix A. If UPLO = 'U', the
! 89: *> leading N-by-N upper triangular part of A contains the
! 90: *> upper triangular part of the matrix A. If UPLO = 'L',
! 91: *> the leading N-by-N lower triangular part of A contains
! 92: *> the lower triangular part of the matrix A.
! 93: *>
! 94: *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
! 95: *> matrix Z of eigenvectors. The eigenvectors are normalized
! 96: *> as follows:
! 97: *> if ITYPE = 1 or 2, Z**T*B*Z = I;
! 98: *> if ITYPE = 3, Z**T*inv(B)*Z = I.
! 99: *> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
! 100: *> or the lower triangle (if UPLO='L') of A, including the
! 101: *> diagonal, is destroyed.
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in] LDA
! 105: *> \verbatim
! 106: *> LDA is INTEGER
! 107: *> The leading dimension of the array A. LDA >= max(1,N).
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[in,out] B
! 111: *> \verbatim
! 112: *> B is DOUBLE PRECISION array, dimension (LDB, N)
! 113: *> On entry, the symmetric matrix B. If UPLO = 'U', the
! 114: *> leading N-by-N upper triangular part of B contains the
! 115: *> upper triangular part of the matrix B. If UPLO = 'L',
! 116: *> the leading N-by-N lower triangular part of B contains
! 117: *> the lower triangular part of the matrix B.
! 118: *>
! 119: *> On exit, if INFO <= N, the part of B containing the matrix is
! 120: *> overwritten by the triangular factor U or L from the Cholesky
! 121: *> factorization B = U**T*U or B = L*L**T.
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[in] LDB
! 125: *> \verbatim
! 126: *> LDB is INTEGER
! 127: *> The leading dimension of the array B. LDB >= max(1,N).
! 128: *> \endverbatim
! 129: *>
! 130: *> \param[out] W
! 131: *> \verbatim
! 132: *> W is DOUBLE PRECISION array, dimension (N)
! 133: *> If INFO = 0, the eigenvalues in ascending order.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[out] WORK
! 137: *> \verbatim
! 138: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 139: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[in] LWORK
! 143: *> \verbatim
! 144: *> LWORK is INTEGER
! 145: *> The dimension of the array WORK.
! 146: *> If N <= 1, LWORK >= 1.
! 147: *> If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
! 148: *> If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
! 149: *>
! 150: *> If LWORK = -1, then a workspace query is assumed; the routine
! 151: *> only calculates the optimal sizes of the WORK and IWORK
! 152: *> arrays, returns these values as the first entries of the WORK
! 153: *> and IWORK arrays, and no error message related to LWORK or
! 154: *> LIWORK is issued by XERBLA.
! 155: *> \endverbatim
! 156: *>
! 157: *> \param[out] IWORK
! 158: *> \verbatim
! 159: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
! 160: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
! 161: *> \endverbatim
! 162: *>
! 163: *> \param[in] LIWORK
! 164: *> \verbatim
! 165: *> LIWORK is INTEGER
! 166: *> The dimension of the array IWORK.
! 167: *> If N <= 1, LIWORK >= 1.
! 168: *> If JOBZ = 'N' and N > 1, LIWORK >= 1.
! 169: *> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
! 170: *>
! 171: *> If LIWORK = -1, then a workspace query is assumed; the
! 172: *> routine only calculates the optimal sizes of the WORK and
! 173: *> IWORK arrays, returns these values as the first entries of
! 174: *> the WORK and IWORK arrays, and no error message related to
! 175: *> LWORK or LIWORK is issued by XERBLA.
! 176: *> \endverbatim
! 177: *>
! 178: *> \param[out] INFO
! 179: *> \verbatim
! 180: *> INFO is INTEGER
! 181: *> = 0: successful exit
! 182: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 183: *> > 0: DPOTRF or DSYEVD returned an error code:
! 184: *> <= N: if INFO = i and JOBZ = 'N', then the algorithm
! 185: *> failed to converge; i off-diagonal elements of an
! 186: *> intermediate tridiagonal form did not converge to
! 187: *> zero;
! 188: *> if INFO = i and JOBZ = 'V', then the algorithm
! 189: *> failed to compute an eigenvalue while working on
! 190: *> the submatrix lying in rows and columns INFO/(N+1)
! 191: *> through mod(INFO,N+1);
! 192: *> > N: if INFO = N + i, for 1 <= i <= N, then the leading
! 193: *> minor of order i of B is not positive definite.
! 194: *> The factorization of B could not be completed and
! 195: *> no eigenvalues or eigenvectors were computed.
! 196: *> \endverbatim
! 197: *
! 198: * Authors:
! 199: * ========
! 200: *
! 201: *> \author Univ. of Tennessee
! 202: *> \author Univ. of California Berkeley
! 203: *> \author Univ. of Colorado Denver
! 204: *> \author NAG Ltd.
! 205: *
! 206: *> \date November 2011
! 207: *
! 208: *> \ingroup doubleSYeigen
! 209: *
! 210: *> \par Further Details:
! 211: * =====================
! 212: *>
! 213: *> \verbatim
! 214: *>
! 215: *> Modified so that no backsubstitution is performed if DSYEVD fails to
! 216: *> converge (NEIG in old code could be greater than N causing out of
! 217: *> bounds reference to A - reported by Ralf Meyer). Also corrected the
! 218: *> description of INFO and the test on ITYPE. Sven, 16 Feb 05.
! 219: *> \endverbatim
! 220: *
! 221: *> \par Contributors:
! 222: * ==================
! 223: *>
! 224: *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
! 225: *>
! 226: * =====================================================================
1.1 bertrand 227: SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
228: $ LWORK, IWORK, LIWORK, INFO )
229: *
1.9 ! bertrand 230: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 231: * -- LAPACK is a software package provided by Univ. of Tennessee, --
232: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 233: * November 2011
1.1 bertrand 234: *
235: * .. Scalar Arguments ..
236: CHARACTER JOBZ, UPLO
237: INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
238: * ..
239: * .. Array Arguments ..
240: INTEGER IWORK( * )
241: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
242: * ..
243: *
244: * =====================================================================
245: *
246: * .. Parameters ..
247: DOUBLE PRECISION ONE
248: PARAMETER ( ONE = 1.0D+0 )
249: * ..
250: * .. Local Scalars ..
251: LOGICAL LQUERY, UPPER, WANTZ
252: CHARACTER TRANS
253: INTEGER LIOPT, LIWMIN, LOPT, LWMIN
254: * ..
255: * .. External Functions ..
256: LOGICAL LSAME
257: EXTERNAL LSAME
258: * ..
259: * .. External Subroutines ..
260: EXTERNAL DPOTRF, DSYEVD, DSYGST, DTRMM, DTRSM, XERBLA
261: * ..
262: * .. Intrinsic Functions ..
263: INTRINSIC DBLE, MAX
264: * ..
265: * .. Executable Statements ..
266: *
267: * Test the input parameters.
268: *
269: WANTZ = LSAME( JOBZ, 'V' )
270: UPPER = LSAME( UPLO, 'U' )
271: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
272: *
273: INFO = 0
274: IF( N.LE.1 ) THEN
275: LIWMIN = 1
276: LWMIN = 1
277: ELSE IF( WANTZ ) THEN
278: LIWMIN = 3 + 5*N
279: LWMIN = 1 + 6*N + 2*N**2
280: ELSE
281: LIWMIN = 1
282: LWMIN = 2*N + 1
283: END IF
284: LOPT = LWMIN
285: LIOPT = LIWMIN
286: IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
287: INFO = -1
288: ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
289: INFO = -2
290: ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
291: INFO = -3
292: ELSE IF( N.LT.0 ) THEN
293: INFO = -4
294: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
295: INFO = -6
296: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
297: INFO = -8
298: END IF
299: *
300: IF( INFO.EQ.0 ) THEN
301: WORK( 1 ) = LOPT
302: IWORK( 1 ) = LIOPT
303: *
304: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
305: INFO = -11
306: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
307: INFO = -13
308: END IF
309: END IF
310: *
311: IF( INFO.NE.0 ) THEN
312: CALL XERBLA( 'DSYGVD', -INFO )
313: RETURN
314: ELSE IF( LQUERY ) THEN
315: RETURN
316: END IF
317: *
318: * Quick return if possible
319: *
320: IF( N.EQ.0 )
321: $ RETURN
322: *
323: * Form a Cholesky factorization of B.
324: *
325: CALL DPOTRF( UPLO, N, B, LDB, INFO )
326: IF( INFO.NE.0 ) THEN
327: INFO = N + INFO
328: RETURN
329: END IF
330: *
331: * Transform problem to standard eigenvalue problem and solve.
332: *
333: CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
334: CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
335: $ INFO )
336: LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
337: LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
338: *
339: IF( WANTZ .AND. INFO.EQ.0 ) THEN
340: *
341: * Backtransform eigenvectors to the original problem.
342: *
343: IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
344: *
345: * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
1.8 bertrand 346: * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
1.1 bertrand 347: *
348: IF( UPPER ) THEN
349: TRANS = 'N'
350: ELSE
351: TRANS = 'T'
352: END IF
353: *
354: CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
355: $ B, LDB, A, LDA )
356: *
357: ELSE IF( ITYPE.EQ.3 ) THEN
358: *
359: * For B*A*x=(lambda)*x;
1.8 bertrand 360: * backtransform eigenvectors: x = L*y or U**T*y
1.1 bertrand 361: *
362: IF( UPPER ) THEN
363: TRANS = 'T'
364: ELSE
365: TRANS = 'N'
366: END IF
367: *
368: CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
369: $ B, LDB, A, LDA )
370: END IF
371: END IF
372: *
373: WORK( 1 ) = LOPT
374: IWORK( 1 ) = LIOPT
375: *
376: RETURN
377: *
378: * End of DSYGVD
379: *
380: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsygvd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack