1: *> \brief \b DSYGVD
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: *> \ingroup doubleSYeigen
207: *
208: *> \par Further Details:
209: * =====================
210: *>
211: *> \verbatim
212: *>
213: *> Modified so that no backsubstitution is performed if DSYEVD fails to
214: *> converge (NEIG in old code could be greater than N causing out of
215: *> bounds reference to A - reported by Ralf Meyer). Also corrected the
216: *> description of INFO and the test on ITYPE. Sven, 16 Feb 05.
217: *> \endverbatim
218: *
219: *> \par Contributors:
220: * ==================
221: *>
222: *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
223: *>
224: * =====================================================================
225: SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
226: $ LWORK, IWORK, LIWORK, INFO )
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 JOBZ, UPLO
234: INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
235: * ..
236: * .. Array Arguments ..
237: INTEGER IWORK( * )
238: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
239: * ..
240: *
241: * =====================================================================
242: *
243: * .. Parameters ..
244: DOUBLE PRECISION ONE
245: PARAMETER ( ONE = 1.0D+0 )
246: * ..
247: * .. Local Scalars ..
248: LOGICAL LQUERY, UPPER, WANTZ
249: CHARACTER TRANS
250: INTEGER LIOPT, LIWMIN, LOPT, LWMIN
251: * ..
252: * .. External Functions ..
253: LOGICAL LSAME
254: EXTERNAL LSAME
255: * ..
256: * .. External Subroutines ..
257: EXTERNAL DPOTRF, DSYEVD, DSYGST, DTRMM, DTRSM, XERBLA
258: * ..
259: * .. Intrinsic Functions ..
260: INTRINSIC DBLE, MAX
261: * ..
262: * .. Executable Statements ..
263: *
264: * Test the input parameters.
265: *
266: WANTZ = LSAME( JOBZ, 'V' )
267: UPPER = LSAME( UPLO, 'U' )
268: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
269: *
270: INFO = 0
271: IF( N.LE.1 ) THEN
272: LIWMIN = 1
273: LWMIN = 1
274: ELSE IF( WANTZ ) THEN
275: LIWMIN = 3 + 5*N
276: LWMIN = 1 + 6*N + 2*N**2
277: ELSE
278: LIWMIN = 1
279: LWMIN = 2*N + 1
280: END IF
281: LOPT = LWMIN
282: LIOPT = LIWMIN
283: IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
284: INFO = -1
285: ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
286: INFO = -2
287: ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
288: INFO = -3
289: ELSE IF( N.LT.0 ) THEN
290: INFO = -4
291: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
292: INFO = -6
293: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
294: INFO = -8
295: END IF
296: *
297: IF( INFO.EQ.0 ) THEN
298: WORK( 1 ) = LOPT
299: IWORK( 1 ) = LIOPT
300: *
301: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
302: INFO = -11
303: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
304: INFO = -13
305: END IF
306: END IF
307: *
308: IF( INFO.NE.0 ) THEN
309: CALL XERBLA( 'DSYGVD', -INFO )
310: RETURN
311: ELSE IF( LQUERY ) THEN
312: RETURN
313: END IF
314: *
315: * Quick return if possible
316: *
317: IF( N.EQ.0 )
318: $ RETURN
319: *
320: * Form a Cholesky factorization of B.
321: *
322: CALL DPOTRF( UPLO, N, B, LDB, INFO )
323: IF( INFO.NE.0 ) THEN
324: INFO = N + INFO
325: RETURN
326: END IF
327: *
328: * Transform problem to standard eigenvalue problem and solve.
329: *
330: CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
331: CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
332: $ INFO )
333: LOPT = INT( MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) )
334: LIOPT = INT( MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) )
335: *
336: IF( WANTZ .AND. INFO.EQ.0 ) THEN
337: *
338: * Backtransform eigenvectors to the original problem.
339: *
340: IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
341: *
342: * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
343: * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
344: *
345: IF( UPPER ) THEN
346: TRANS = 'N'
347: ELSE
348: TRANS = 'T'
349: END IF
350: *
351: CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
352: $ B, LDB, A, LDA )
353: *
354: ELSE IF( ITYPE.EQ.3 ) THEN
355: *
356: * For B*A*x=(lambda)*x;
357: * backtransform eigenvectors: x = L*y or U**T*y
358: *
359: IF( UPPER ) THEN
360: TRANS = 'T'
361: ELSE
362: TRANS = 'N'
363: END IF
364: *
365: CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE,
366: $ B, LDB, A, LDA )
367: END IF
368: END IF
369: *
370: WORK( 1 ) = LOPT
371: IWORK( 1 ) = LIOPT
372: *
373: RETURN
374: *
375: * End of DSYGVD
376: *
377: 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