1: *> \brief \b DSYGVX
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYGVX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygvx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygvx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygvx.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
22: * VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK,
23: * LWORK, IWORK, IFAIL, INFO )
24: *
25: * .. Scalar Arguments ..
26: * CHARACTER JOBZ, RANGE, UPLO
27: * INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
28: * DOUBLE PRECISION ABSTOL, VL, VU
29: * ..
30: * .. Array Arguments ..
31: * INTEGER IFAIL( * ), IWORK( * )
32: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
33: * $ Z( LDZ, * )
34: * ..
35: *
36: *
37: *> \par Purpose:
38: * =============
39: *>
40: *> \verbatim
41: *>
42: *> DSYGVX computes selected eigenvalues, and optionally, eigenvectors
43: *> of a real generalized symmetric-definite eigenproblem, of the form
44: *> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A
45: *> and B are assumed to be symmetric and B is also positive definite.
46: *> Eigenvalues and eigenvectors can be selected by specifying either a
47: *> range of values or a range of indices for the desired eigenvalues.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] ITYPE
54: *> \verbatim
55: *> ITYPE is INTEGER
56: *> Specifies the problem type to be solved:
57: *> = 1: A*x = (lambda)*B*x
58: *> = 2: A*B*x = (lambda)*x
59: *> = 3: B*A*x = (lambda)*x
60: *> \endverbatim
61: *>
62: *> \param[in] JOBZ
63: *> \verbatim
64: *> JOBZ is CHARACTER*1
65: *> = 'N': Compute eigenvalues only;
66: *> = 'V': Compute eigenvalues and eigenvectors.
67: *> \endverbatim
68: *>
69: *> \param[in] RANGE
70: *> \verbatim
71: *> RANGE is CHARACTER*1
72: *> = 'A': all eigenvalues will be found.
73: *> = 'V': all eigenvalues in the half-open interval (VL,VU]
74: *> will be found.
75: *> = 'I': the IL-th through IU-th eigenvalues will be found.
76: *> \endverbatim
77: *>
78: *> \param[in] UPLO
79: *> \verbatim
80: *> UPLO is CHARACTER*1
81: *> = 'U': Upper triangle of A and B are stored;
82: *> = 'L': Lower triangle of A and B are stored.
83: *> \endverbatim
84: *>
85: *> \param[in] N
86: *> \verbatim
87: *> N is INTEGER
88: *> The order of the matrix pencil (A,B). N >= 0.
89: *> \endverbatim
90: *>
91: *> \param[in,out] A
92: *> \verbatim
93: *> A is DOUBLE PRECISION array, dimension (LDA, N)
94: *> On entry, the symmetric matrix A. If UPLO = 'U', the
95: *> leading N-by-N upper triangular part of A contains the
96: *> upper triangular part of the matrix A. If UPLO = 'L',
97: *> the leading N-by-N lower triangular part of A contains
98: *> the lower triangular part of the matrix A.
99: *>
100: *> On exit, the lower triangle (if UPLO='L') or the upper
101: *> triangle (if UPLO='U') of A, including the diagonal, is
102: *> destroyed.
103: *> \endverbatim
104: *>
105: *> \param[in] LDA
106: *> \verbatim
107: *> LDA is INTEGER
108: *> The leading dimension of the array A. LDA >= max(1,N).
109: *> \endverbatim
110: *>
111: *> \param[in,out] B
112: *> \verbatim
113: *> B is DOUBLE PRECISION array, dimension (LDB, N)
114: *> On entry, the symmetric matrix B. If UPLO = 'U', the
115: *> leading N-by-N upper triangular part of B contains the
116: *> upper triangular part of the matrix B. If UPLO = 'L',
117: *> the leading N-by-N lower triangular part of B contains
118: *> the lower triangular part of the matrix B.
119: *>
120: *> On exit, if INFO <= N, the part of B containing the matrix is
121: *> overwritten by the triangular factor U or L from the Cholesky
122: *> factorization B = U**T*U or B = L*L**T.
123: *> \endverbatim
124: *>
125: *> \param[in] LDB
126: *> \verbatim
127: *> LDB is INTEGER
128: *> The leading dimension of the array B. LDB >= max(1,N).
129: *> \endverbatim
130: *>
131: *> \param[in] VL
132: *> \verbatim
133: *> VL is DOUBLE PRECISION
134: *> If RANGE='V', the lower bound of the interval to
135: *> be searched for eigenvalues. VL < VU.
136: *> Not referenced if RANGE = 'A' or 'I'.
137: *> \endverbatim
138: *>
139: *> \param[in] VU
140: *> \verbatim
141: *> VU is DOUBLE PRECISION
142: *> If RANGE='V', the upper bound of the interval to
143: *> be searched for eigenvalues. VL < VU.
144: *> Not referenced if RANGE = 'A' or 'I'.
145: *> \endverbatim
146: *>
147: *> \param[in] IL
148: *> \verbatim
149: *> IL is INTEGER
150: *> If RANGE='I', the index of the
151: *> smallest eigenvalue to be returned.
152: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
153: *> Not referenced if RANGE = 'A' or 'V'.
154: *> \endverbatim
155: *>
156: *> \param[in] IU
157: *> \verbatim
158: *> IU is INTEGER
159: *> If RANGE='I', the index of the
160: *> largest eigenvalue to be returned.
161: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
162: *> Not referenced if RANGE = 'A' or 'V'.
163: *> \endverbatim
164: *>
165: *> \param[in] ABSTOL
166: *> \verbatim
167: *> ABSTOL is DOUBLE PRECISION
168: *> The absolute error tolerance for the eigenvalues.
169: *> An approximate eigenvalue is accepted as converged
170: *> when it is determined to lie in an interval [a,b]
171: *> of width less than or equal to
172: *>
173: *> ABSTOL + EPS * max( |a|,|b| ) ,
174: *>
175: *> where EPS is the machine precision. If ABSTOL is less than
176: *> or equal to zero, then EPS*|T| will be used in its place,
177: *> where |T| is the 1-norm of the tridiagonal matrix obtained
178: *> by reducing C to tridiagonal form, where C is the symmetric
179: *> matrix of the standard symmetric problem to which the
180: *> generalized problem is transformed.
181: *>
182: *> Eigenvalues will be computed most accurately when ABSTOL is
183: *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
184: *> If this routine returns with INFO>0, indicating that some
185: *> eigenvectors did not converge, try setting ABSTOL to
186: *> 2*DLAMCH('S').
187: *> \endverbatim
188: *>
189: *> \param[out] M
190: *> \verbatim
191: *> M is INTEGER
192: *> The total number of eigenvalues found. 0 <= M <= N.
193: *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
194: *> \endverbatim
195: *>
196: *> \param[out] W
197: *> \verbatim
198: *> W is DOUBLE PRECISION array, dimension (N)
199: *> On normal exit, the first M elements contain the selected
200: *> eigenvalues in ascending order.
201: *> \endverbatim
202: *>
203: *> \param[out] Z
204: *> \verbatim
205: *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
206: *> If JOBZ = 'N', then Z is not referenced.
207: *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
208: *> contain the orthonormal eigenvectors of the matrix A
209: *> corresponding to the selected eigenvalues, with the i-th
210: *> column of Z holding the eigenvector associated with W(i).
211: *> The eigenvectors are normalized as follows:
212: *> if ITYPE = 1 or 2, Z**T*B*Z = I;
213: *> if ITYPE = 3, Z**T*inv(B)*Z = I.
214: *>
215: *> If an eigenvector fails to converge, then that column of Z
216: *> contains the latest approximation to the eigenvector, and the
217: *> index of the eigenvector is returned in IFAIL.
218: *> Note: the user must ensure that at least max(1,M) columns are
219: *> supplied in the array Z; if RANGE = 'V', the exact value of M
220: *> is not known in advance and an upper bound must be used.
221: *> \endverbatim
222: *>
223: *> \param[in] LDZ
224: *> \verbatim
225: *> LDZ is INTEGER
226: *> The leading dimension of the array Z. LDZ >= 1, and if
227: *> JOBZ = 'V', LDZ >= max(1,N).
228: *> \endverbatim
229: *>
230: *> \param[out] WORK
231: *> \verbatim
232: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
233: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
234: *> \endverbatim
235: *>
236: *> \param[in] LWORK
237: *> \verbatim
238: *> LWORK is INTEGER
239: *> The length of the array WORK. LWORK >= max(1,8*N).
240: *> For optimal efficiency, LWORK >= (NB+3)*N,
241: *> where NB is the blocksize for DSYTRD returned by ILAENV.
242: *>
243: *> If LWORK = -1, then a workspace query is assumed; the routine
244: *> only calculates the optimal size of the WORK array, returns
245: *> this value as the first entry of the WORK array, and no error
246: *> message related to LWORK is issued by XERBLA.
247: *> \endverbatim
248: *>
249: *> \param[out] IWORK
250: *> \verbatim
251: *> IWORK is INTEGER array, dimension (5*N)
252: *> \endverbatim
253: *>
254: *> \param[out] IFAIL
255: *> \verbatim
256: *> IFAIL is INTEGER array, dimension (N)
257: *> If JOBZ = 'V', then if INFO = 0, the first M elements of
258: *> IFAIL are zero. If INFO > 0, then IFAIL contains the
259: *> indices of the eigenvectors that failed to converge.
260: *> If JOBZ = 'N', then IFAIL is not referenced.
261: *> \endverbatim
262: *>
263: *> \param[out] INFO
264: *> \verbatim
265: *> INFO is INTEGER
266: *> = 0: successful exit
267: *> < 0: if INFO = -i, the i-th argument had an illegal value
268: *> > 0: DPOTRF or DSYEVX returned an error code:
269: *> <= N: if INFO = i, DSYEVX failed to converge;
270: *> i eigenvectors failed to converge. Their indices
271: *> are stored in array IFAIL.
272: *> > N: if INFO = N + i, for 1 <= i <= N, then the leading
273: *> minor of order i of B is not positive definite.
274: *> The factorization of B could not be completed and
275: *> no eigenvalues or eigenvectors were computed.
276: *> \endverbatim
277: *
278: * Authors:
279: * ========
280: *
281: *> \author Univ. of Tennessee
282: *> \author Univ. of California Berkeley
283: *> \author Univ. of Colorado Denver
284: *> \author NAG Ltd.
285: *
286: *> \ingroup doubleSYeigen
287: *
288: *> \par Contributors:
289: * ==================
290: *>
291: *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
292: *
293: * =====================================================================
294: SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
295: $ VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK,
296: $ LWORK, IWORK, IFAIL, INFO )
297: *
298: * -- LAPACK driver routine --
299: * -- LAPACK is a software package provided by Univ. of Tennessee, --
300: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301: *
302: * .. Scalar Arguments ..
303: CHARACTER JOBZ, RANGE, UPLO
304: INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
305: DOUBLE PRECISION ABSTOL, VL, VU
306: * ..
307: * .. Array Arguments ..
308: INTEGER IFAIL( * ), IWORK( * )
309: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
310: $ Z( LDZ, * )
311: * ..
312: *
313: * =====================================================================
314: *
315: * .. Parameters ..
316: DOUBLE PRECISION ONE
317: PARAMETER ( ONE = 1.0D+0 )
318: * ..
319: * .. Local Scalars ..
320: LOGICAL ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ
321: CHARACTER TRANS
322: INTEGER LWKMIN, LWKOPT, NB
323: * ..
324: * .. External Functions ..
325: LOGICAL LSAME
326: INTEGER ILAENV
327: EXTERNAL LSAME, ILAENV
328: * ..
329: * .. External Subroutines ..
330: EXTERNAL DPOTRF, DSYEVX, DSYGST, DTRMM, DTRSM, XERBLA
331: * ..
332: * .. Intrinsic Functions ..
333: INTRINSIC MAX, MIN
334: * ..
335: * .. Executable Statements ..
336: *
337: * Test the input parameters.
338: *
339: UPPER = LSAME( UPLO, 'U' )
340: WANTZ = LSAME( JOBZ, 'V' )
341: ALLEIG = LSAME( RANGE, 'A' )
342: VALEIG = LSAME( RANGE, 'V' )
343: INDEIG = LSAME( RANGE, 'I' )
344: LQUERY = ( LWORK.EQ.-1 )
345: *
346: INFO = 0
347: IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
348: INFO = -1
349: ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
350: INFO = -2
351: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
352: INFO = -3
353: ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
354: INFO = -4
355: ELSE IF( N.LT.0 ) THEN
356: INFO = -5
357: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
358: INFO = -7
359: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
360: INFO = -9
361: ELSE
362: IF( VALEIG ) THEN
363: IF( N.GT.0 .AND. VU.LE.VL )
364: $ INFO = -11
365: ELSE IF( INDEIG ) THEN
366: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
367: INFO = -12
368: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
369: INFO = -13
370: END IF
371: END IF
372: END IF
373: IF (INFO.EQ.0) THEN
374: IF (LDZ.LT.1 .OR. (WANTZ .AND. LDZ.LT.N)) THEN
375: INFO = -18
376: END IF
377: END IF
378: *
379: IF( INFO.EQ.0 ) THEN
380: LWKMIN = MAX( 1, 8*N )
381: NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
382: LWKOPT = MAX( LWKMIN, ( NB + 3 )*N )
383: WORK( 1 ) = LWKOPT
384: *
385: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
386: INFO = -20
387: END IF
388: END IF
389: *
390: IF( INFO.NE.0 ) THEN
391: CALL XERBLA( 'DSYGVX', -INFO )
392: RETURN
393: ELSE IF( LQUERY ) THEN
394: RETURN
395: END IF
396: *
397: * Quick return if possible
398: *
399: M = 0
400: IF( N.EQ.0 ) THEN
401: RETURN
402: END IF
403: *
404: * Form a Cholesky factorization of B.
405: *
406: CALL DPOTRF( UPLO, N, B, LDB, INFO )
407: IF( INFO.NE.0 ) THEN
408: INFO = N + INFO
409: RETURN
410: END IF
411: *
412: * Transform problem to standard eigenvalue problem and solve.
413: *
414: CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
415: CALL DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL,
416: $ M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
417: *
418: IF( WANTZ ) THEN
419: *
420: * Backtransform eigenvectors to the original problem.
421: *
422: IF( INFO.GT.0 )
423: $ M = INFO - 1
424: IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
425: *
426: * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
427: * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
428: *
429: IF( UPPER ) THEN
430: TRANS = 'N'
431: ELSE
432: TRANS = 'T'
433: END IF
434: *
435: CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
436: $ LDB, Z, LDZ )
437: *
438: ELSE IF( ITYPE.EQ.3 ) THEN
439: *
440: * For B*A*x=(lambda)*x;
441: * backtransform eigenvectors: x = L*y or U**T*y
442: *
443: IF( UPPER ) THEN
444: TRANS = 'T'
445: ELSE
446: TRANS = 'N'
447: END IF
448: *
449: CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
450: $ LDB, Z, LDZ )
451: END IF
452: END IF
453: *
454: * Set WORK(1) to optimal workspace size.
455: *
456: WORK( 1 ) = LWKOPT
457: *
458: RETURN
459: *
460: * End of DSYGVX
461: *
462: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsygvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack