LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
dsygvx.f
Go to the documentation of this file.
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 *> \endverbatim
135 *>
136 *> \param[in] VU
137 *> \verbatim
138 *> VU is DOUBLE PRECISION
139 *> If RANGE='V', the lower and upper bounds of the interval to
140 *> be searched for eigenvalues. VL < VU.
141 *> Not referenced if RANGE = 'A' or 'I'.
142 *> \endverbatim
143 *>
144 *> \param[in] IL
145 *> \verbatim
146 *> IL is INTEGER
147 *> \endverbatim
148 *>
149 *> \param[in] IU
150 *> \verbatim
151 *> IU is INTEGER
152 *> If RANGE='I', the indices (in ascending order) of the
153 *> smallest and largest eigenvalues to be returned.
154 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
155 *> Not referenced if RANGE = 'A' or 'V'.
156 *> \endverbatim
157 *>
158 *> \param[in] ABSTOL
159 *> \verbatim
160 *> ABSTOL is DOUBLE PRECISION
161 *> The absolute error tolerance for the eigenvalues.
162 *> An approximate eigenvalue is accepted as converged
163 *> when it is determined to lie in an interval [a,b]
164 *> of width less than or equal to
165 *>
166 *> ABSTOL + EPS * max( |a|,|b| ) ,
167 *>
168 *> where EPS is the machine precision. If ABSTOL is less than
169 *> or equal to zero, then EPS*|T| will be used in its place,
170 *> where |T| is the 1-norm of the tridiagonal matrix obtained
171 *> by reducing C to tridiagonal form, where C is the symmetric
172 *> matrix of the standard symmetric problem to which the
173 *> generalized problem is transformed.
174 *>
175 *> Eigenvalues will be computed most accurately when ABSTOL is
176 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
177 *> If this routine returns with INFO>0, indicating that some
178 *> eigenvectors did not converge, try setting ABSTOL to
179 *> 2*DLAMCH('S').
180 *> \endverbatim
181 *>
182 *> \param[out] M
183 *> \verbatim
184 *> M is INTEGER
185 *> The total number of eigenvalues found. 0 <= M <= N.
186 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
187 *> \endverbatim
188 *>
189 *> \param[out] W
190 *> \verbatim
191 *> W is DOUBLE PRECISION array, dimension (N)
192 *> On normal exit, the first M elements contain the selected
193 *> eigenvalues in ascending order.
194 *> \endverbatim
195 *>
196 *> \param[out] Z
197 *> \verbatim
198 *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
199 *> If JOBZ = 'N', then Z is not referenced.
200 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
201 *> contain the orthonormal eigenvectors of the matrix A
202 *> corresponding to the selected eigenvalues, with the i-th
203 *> column of Z holding the eigenvector associated with W(i).
204 *> The eigenvectors are normalized as follows:
205 *> if ITYPE = 1 or 2, Z**T*B*Z = I;
206 *> if ITYPE = 3, Z**T*inv(B)*Z = I.
207 *>
208 *> If an eigenvector fails to converge, then that column of Z
209 *> contains the latest approximation to the eigenvector, and the
210 *> index of the eigenvector is returned in IFAIL.
211 *> Note: the user must ensure that at least max(1,M) columns are
212 *> supplied in the array Z; if RANGE = 'V', the exact value of M
213 *> is not known in advance and an upper bound must be used.
214 *> \endverbatim
215 *>
216 *> \param[in] LDZ
217 *> \verbatim
218 *> LDZ is INTEGER
219 *> The leading dimension of the array Z. LDZ >= 1, and if
220 *> JOBZ = 'V', LDZ >= max(1,N).
221 *> \endverbatim
222 *>
223 *> \param[out] WORK
224 *> \verbatim
225 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
226 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
227 *> \endverbatim
228 *>
229 *> \param[in] LWORK
230 *> \verbatim
231 *> LWORK is INTEGER
232 *> The length of the array WORK. LWORK >= max(1,8*N).
233 *> For optimal efficiency, LWORK >= (NB+3)*N,
234 *> where NB is the blocksize for DSYTRD returned by ILAENV.
235 *>
236 *> If LWORK = -1, then a workspace query is assumed; the routine
237 *> only calculates the optimal size of the WORK array, returns
238 *> this value as the first entry of the WORK array, and no error
239 *> message related to LWORK is issued by XERBLA.
240 *> \endverbatim
241 *>
242 *> \param[out] IWORK
243 *> \verbatim
244 *> IWORK is INTEGER array, dimension (5*N)
245 *> \endverbatim
246 *>
247 *> \param[out] IFAIL
248 *> \verbatim
249 *> IFAIL is INTEGER array, dimension (N)
250 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
251 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
252 *> indices of the eigenvectors that failed to converge.
253 *> If JOBZ = 'N', then IFAIL is not referenced.
254 *> \endverbatim
255 *>
256 *> \param[out] INFO
257 *> \verbatim
258 *> INFO is INTEGER
259 *> = 0: successful exit
260 *> < 0: if INFO = -i, the i-th argument had an illegal value
261 *> > 0: DPOTRF or DSYEVX returned an error code:
262 *> <= N: if INFO = i, DSYEVX failed to converge;
263 *> i eigenvectors failed to converge. Their indices
264 *> are stored in array IFAIL.
265 *> > N: if INFO = N + i, for 1 <= i <= N, then the leading
266 *> minor of order i of B is not positive definite.
267 *> The factorization of B could not be completed and
268 *> no eigenvalues or eigenvectors were computed.
269 *> \endverbatim
270 *
271 * Authors:
272 * ========
273 *
274 *> \author Univ. of Tennessee
275 *> \author Univ. of California Berkeley
276 *> \author Univ. of Colorado Denver
277 *> \author NAG Ltd.
278 *
279 *> \date November 2015
280 *
281 *> \ingroup doubleSYeigen
282 *
283 *> \par Contributors:
284 * ==================
285 *>
286 *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
287 *
288 * =====================================================================
289  SUBROUTINE dsygvx( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
290  $ vl, vu, il, iu, abstol, m, w, z, ldz, work,
291  $ lwork, iwork, ifail, info )
292 *
293 * -- LAPACK driver routine (version 3.6.0) --
294 * -- LAPACK is a software package provided by Univ. of Tennessee, --
295 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
296 * November 2015
297 *
298 * .. Scalar Arguments ..
299  CHARACTER JOBZ, RANGE, UPLO
300  INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
301  DOUBLE PRECISION ABSTOL, VL, VU
302 * ..
303 * .. Array Arguments ..
304  INTEGER IFAIL( * ), IWORK( * )
305  DOUBLE PRECISION A( lda, * ), B( ldb, * ), W( * ), WORK( * ),
306  $ z( ldz, * )
307 * ..
308 *
309 * =====================================================================
310 *
311 * .. Parameters ..
312  DOUBLE PRECISION ONE
313  parameter( one = 1.0d+0 )
314 * ..
315 * .. Local Scalars ..
316  LOGICAL ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ
317  CHARACTER TRANS
318  INTEGER LWKMIN, LWKOPT, NB
319 * ..
320 * .. External Functions ..
321  LOGICAL LSAME
322  INTEGER ILAENV
323  EXTERNAL lsame, ilaenv
324 * ..
325 * .. External Subroutines ..
326  EXTERNAL dpotrf, dsyevx, dsygst, dtrmm, dtrsm, xerbla
327 * ..
328 * .. Intrinsic Functions ..
329  INTRINSIC max, min
330 * ..
331 * .. Executable Statements ..
332 *
333 * Test the input parameters.
334 *
335  upper = lsame( uplo, 'U' )
336  wantz = lsame( jobz, 'V' )
337  alleig = lsame( range, 'A' )
338  valeig = lsame( range, 'V' )
339  indeig = lsame( range, 'I' )
340  lquery = ( lwork.EQ.-1 )
341 *
342  info = 0
343  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
344  info = -1
345  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
346  info = -2
347  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
348  info = -3
349  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
350  info = -4
351  ELSE IF( n.LT.0 ) THEN
352  info = -5
353  ELSE IF( lda.LT.max( 1, n ) ) THEN
354  info = -7
355  ELSE IF( ldb.LT.max( 1, n ) ) THEN
356  info = -9
357  ELSE
358  IF( valeig ) THEN
359  IF( n.GT.0 .AND. vu.LE.vl )
360  $ info = -11
361  ELSE IF( indeig ) THEN
362  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
363  info = -12
364  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
365  info = -13
366  END IF
367  END IF
368  END IF
369  IF (info.EQ.0) THEN
370  IF (ldz.LT.1 .OR. (wantz .AND. ldz.LT.n)) THEN
371  info = -18
372  END IF
373  END IF
374 *
375  IF( info.EQ.0 ) THEN
376  lwkmin = max( 1, 8*n )
377  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
378  lwkopt = max( lwkmin, ( nb + 3 )*n )
379  work( 1 ) = lwkopt
380 *
381  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
382  info = -20
383  END IF
384  END IF
385 *
386  IF( info.NE.0 ) THEN
387  CALL xerbla( 'DSYGVX', -info )
388  RETURN
389  ELSE IF( lquery ) THEN
390  RETURN
391  END IF
392 *
393 * Quick return if possible
394 *
395  m = 0
396  IF( n.EQ.0 ) THEN
397  RETURN
398  END IF
399 *
400 * Form a Cholesky factorization of B.
401 *
402  CALL dpotrf( uplo, n, b, ldb, info )
403  IF( info.NE.0 ) THEN
404  info = n + info
405  RETURN
406  END IF
407 *
408 * Transform problem to standard eigenvalue problem and solve.
409 *
410  CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
411  CALL dsyevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,
412  $ m, w, z, ldz, work, lwork, iwork, ifail, info )
413 *
414  IF( wantz ) THEN
415 *
416 * Backtransform eigenvectors to the original problem.
417 *
418  IF( info.GT.0 )
419  $ m = info - 1
420  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
421 *
422 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
423 * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
424 *
425  IF( upper ) THEN
426  trans = 'N'
427  ELSE
428  trans = 'T'
429  END IF
430 *
431  CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, m, one, b,
432  $ ldb, z, ldz )
433 *
434  ELSE IF( itype.EQ.3 ) THEN
435 *
436 * For B*A*x=(lambda)*x;
437 * backtransform eigenvectors: x = L*y or U**T*y
438 *
439  IF( upper ) THEN
440  trans = 'T'
441  ELSE
442  trans = 'N'
443  END IF
444 *
445  CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, m, one, b,
446  $ ldb, z, ldz )
447  END IF
448  END IF
449 *
450 * Set WORK(1) to optimal workspace size.
451 *
452  work( 1 ) = lwkopt
453 *
454  RETURN
455 *
456 * End of DSYGVX
457 *
458  END
subroutine dsyevx(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO)
DSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices ...
Definition: dsyevx.f:248
subroutine dsygst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
DSYGST
Definition: dsygst.f:129
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:183
subroutine dpotrf(UPLO, N, A, LDA, INFO)
DPOTRF
Definition: dpotrf.f:109
subroutine dsygvx(ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO)
DSYGVX
Definition: dsygvx.f:292