Annotation of rpl/lapack/lapack/dsygvx.f, revision 1.19
1.15 bertrand 1: *> \brief \b DSYGVX
1.10 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.18 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.10 bertrand 7: *
8: *> \htmlonly
1.18 bertrand 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">
1.10 bertrand 15: *> [TXT]</a>
1.18 bertrand 16: *> \endhtmlonly
1.10 bertrand 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 )
1.18 bertrand 24: *
1.10 bertrand 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: * ..
1.18 bertrand 35: *
1.10 bertrand 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
1.16 bertrand 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'.
1.10 bertrand 137: *> \endverbatim
138: *>
139: *> \param[in] VU
140: *> \verbatim
141: *> VU is DOUBLE PRECISION
1.16 bertrand 142: *> If RANGE='V', the upper bound of the interval to
1.10 bertrand 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
1.16 bertrand 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'.
1.10 bertrand 154: *> \endverbatim
155: *>
156: *> \param[in] IU
157: *> \verbatim
158: *> IU is INTEGER
1.16 bertrand 159: *> If RANGE='I', the index of the
160: *> largest eigenvalue to be returned.
1.10 bertrand 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: *
1.18 bertrand 281: *> \author Univ. of Tennessee
282: *> \author Univ. of California Berkeley
283: *> \author Univ. of Colorado Denver
284: *> \author NAG Ltd.
1.10 bertrand 285: *
1.16 bertrand 286: *> \date June 2016
1.10 bertrand 287: *
288: *> \ingroup doubleSYeigen
289: *
290: *> \par Contributors:
291: * ==================
292: *>
293: *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
294: *
295: * =====================================================================
1.1 bertrand 296: SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
297: $ VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK,
298: $ LWORK, IWORK, IFAIL, INFO )
299: *
1.18 bertrand 300: * -- LAPACK driver routine (version 3.7.0) --
1.1 bertrand 301: * -- LAPACK is a software package provided by Univ. of Tennessee, --
302: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16 bertrand 303: * June 2016
1.1 bertrand 304: *
305: * .. Scalar Arguments ..
306: CHARACTER JOBZ, RANGE, UPLO
307: INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
308: DOUBLE PRECISION ABSTOL, VL, VU
309: * ..
310: * .. Array Arguments ..
311: INTEGER IFAIL( * ), IWORK( * )
312: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
313: $ Z( LDZ, * )
314: * ..
315: *
316: * =====================================================================
317: *
318: * .. Parameters ..
319: DOUBLE PRECISION ONE
320: PARAMETER ( ONE = 1.0D+0 )
321: * ..
322: * .. Local Scalars ..
323: LOGICAL ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ
324: CHARACTER TRANS
325: INTEGER LWKMIN, LWKOPT, NB
326: * ..
327: * .. External Functions ..
328: LOGICAL LSAME
329: INTEGER ILAENV
330: EXTERNAL LSAME, ILAENV
331: * ..
332: * .. External Subroutines ..
333: EXTERNAL DPOTRF, DSYEVX, DSYGST, DTRMM, DTRSM, XERBLA
334: * ..
335: * .. Intrinsic Functions ..
336: INTRINSIC MAX, MIN
337: * ..
338: * .. Executable Statements ..
339: *
340: * Test the input parameters.
341: *
342: UPPER = LSAME( UPLO, 'U' )
343: WANTZ = LSAME( JOBZ, 'V' )
344: ALLEIG = LSAME( RANGE, 'A' )
345: VALEIG = LSAME( RANGE, 'V' )
346: INDEIG = LSAME( RANGE, 'I' )
347: LQUERY = ( LWORK.EQ.-1 )
348: *
349: INFO = 0
350: IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
351: INFO = -1
352: ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
353: INFO = -2
354: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
355: INFO = -3
356: ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
357: INFO = -4
358: ELSE IF( N.LT.0 ) THEN
359: INFO = -5
360: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
361: INFO = -7
362: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
363: INFO = -9
364: ELSE
365: IF( VALEIG ) THEN
366: IF( N.GT.0 .AND. VU.LE.VL )
367: $ INFO = -11
368: ELSE IF( INDEIG ) THEN
369: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
370: INFO = -12
371: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
372: INFO = -13
373: END IF
374: END IF
375: END IF
376: IF (INFO.EQ.0) THEN
377: IF (LDZ.LT.1 .OR. (WANTZ .AND. LDZ.LT.N)) THEN
378: INFO = -18
379: END IF
380: END IF
381: *
382: IF( INFO.EQ.0 ) THEN
383: LWKMIN = MAX( 1, 8*N )
384: NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
385: LWKOPT = MAX( LWKMIN, ( NB + 3 )*N )
386: WORK( 1 ) = LWKOPT
387: *
388: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
389: INFO = -20
390: END IF
391: END IF
392: *
393: IF( INFO.NE.0 ) THEN
394: CALL XERBLA( 'DSYGVX', -INFO )
395: RETURN
396: ELSE IF( LQUERY ) THEN
397: RETURN
398: END IF
399: *
400: * Quick return if possible
401: *
402: M = 0
403: IF( N.EQ.0 ) THEN
404: RETURN
405: END IF
406: *
407: * Form a Cholesky factorization of B.
408: *
409: CALL DPOTRF( UPLO, N, B, LDB, INFO )
410: IF( INFO.NE.0 ) THEN
411: INFO = N + INFO
412: RETURN
413: END IF
414: *
415: * Transform problem to standard eigenvalue problem and solve.
416: *
417: CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
418: CALL DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL,
419: $ M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
420: *
421: IF( WANTZ ) THEN
422: *
423: * Backtransform eigenvectors to the original problem.
424: *
425: IF( INFO.GT.0 )
426: $ M = INFO - 1
427: IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
428: *
429: * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
1.9 bertrand 430: * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
1.1 bertrand 431: *
432: IF( UPPER ) THEN
433: TRANS = 'N'
434: ELSE
435: TRANS = 'T'
436: END IF
437: *
438: CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
439: $ LDB, Z, LDZ )
440: *
441: ELSE IF( ITYPE.EQ.3 ) THEN
442: *
443: * For B*A*x=(lambda)*x;
1.9 bertrand 444: * backtransform eigenvectors: x = L*y or U**T*y
1.1 bertrand 445: *
446: IF( UPPER ) THEN
447: TRANS = 'T'
448: ELSE
449: TRANS = 'N'
450: END IF
451: *
452: CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
453: $ LDB, Z, LDZ )
454: END IF
455: END IF
456: *
457: * Set WORK(1) to optimal workspace size.
458: *
459: WORK( 1 ) = LWKOPT
460: *
461: RETURN
462: *
463: * End of DSYGVX
464: *
465: 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