LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsygvd.f
Go to the documentation of this file.
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*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] ITYPE
51*> \verbatim
52*> ITYPE is INTEGER
53*> Specifies the problem type to be solved:
54*> = 1: A*x = (lambda)*B*x
55*> = 2: A*B*x = (lambda)*x
56*> = 3: B*A*x = (lambda)*x
57*> \endverbatim
58*>
59*> \param[in] JOBZ
60*> \verbatim
61*> JOBZ is CHARACTER*1
62*> = 'N': Compute eigenvalues only;
63*> = 'V': Compute eigenvalues and eigenvectors.
64*> \endverbatim
65*>
66*> \param[in] UPLO
67*> \verbatim
68*> UPLO is CHARACTER*1
69*> = 'U': Upper triangles of A and B are stored;
70*> = 'L': Lower triangles of A and B are stored.
71*> \endverbatim
72*>
73*> \param[in] N
74*> \verbatim
75*> N is INTEGER
76*> The order of the matrices A and B. N >= 0.
77*> \endverbatim
78*>
79*> \param[in,out] A
80*> \verbatim
81*> A is DOUBLE PRECISION array, dimension (LDA, N)
82*> On entry, the symmetric matrix A. If UPLO = 'U', the
83*> leading N-by-N upper triangular part of A contains the
84*> upper triangular part of the matrix A. If UPLO = 'L',
85*> the leading N-by-N lower triangular part of A contains
86*> the lower triangular part of the matrix A.
87*>
88*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
89*> matrix Z of eigenvectors. The eigenvectors are normalized
90*> as follows:
91*> if ITYPE = 1 or 2, Z**T*B*Z = I;
92*> if ITYPE = 3, Z**T*inv(B)*Z = I.
93*> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
94*> or the lower triangle (if UPLO='L') of A, including the
95*> diagonal, is destroyed.
96*> \endverbatim
97*>
98*> \param[in] LDA
99*> \verbatim
100*> LDA is INTEGER
101*> The leading dimension of the array A. LDA >= max(1,N).
102*> \endverbatim
103*>
104*> \param[in,out] B
105*> \verbatim
106*> B is DOUBLE PRECISION array, dimension (LDB, N)
107*> On entry, the symmetric matrix B. If UPLO = 'U', the
108*> leading N-by-N upper triangular part of B contains the
109*> upper triangular part of the matrix B. If UPLO = 'L',
110*> the leading N-by-N lower triangular part of B contains
111*> the lower triangular part of the matrix B.
112*>
113*> On exit, if INFO <= N, the part of B containing the matrix is
114*> overwritten by the triangular factor U or L from the Cholesky
115*> factorization B = U**T*U or B = L*L**T.
116*> \endverbatim
117*>
118*> \param[in] LDB
119*> \verbatim
120*> LDB is INTEGER
121*> The leading dimension of the array B. LDB >= max(1,N).
122*> \endverbatim
123*>
124*> \param[out] W
125*> \verbatim
126*> W is DOUBLE PRECISION array, dimension (N)
127*> If INFO = 0, the eigenvalues in ascending order.
128*> \endverbatim
129*>
130*> \param[out] WORK
131*> \verbatim
132*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
133*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
134*> \endverbatim
135*>
136*> \param[in] LWORK
137*> \verbatim
138*> LWORK is INTEGER
139*> The dimension of the array WORK.
140*> If N <= 1, LWORK >= 1.
141*> If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
142*> If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
143*>
144*> If LWORK = -1, then a workspace query is assumed; the routine
145*> only calculates the optimal sizes of the WORK and IWORK
146*> arrays, returns these values as the first entries of the WORK
147*> and IWORK arrays, and no error message related to LWORK or
148*> LIWORK is issued by XERBLA.
149*> \endverbatim
150*>
151*> \param[out] IWORK
152*> \verbatim
153*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
154*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
155*> \endverbatim
156*>
157*> \param[in] LIWORK
158*> \verbatim
159*> LIWORK is INTEGER
160*> The dimension of the array IWORK.
161*> If N <= 1, LIWORK >= 1.
162*> If JOBZ = 'N' and N > 1, LIWORK >= 1.
163*> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
164*>
165*> If LIWORK = -1, then a workspace query is assumed; the
166*> routine only calculates the optimal sizes of the WORK and
167*> IWORK arrays, returns these values as the first entries of
168*> the WORK and IWORK arrays, and no error message related to
169*> LWORK or LIWORK is issued by XERBLA.
170*> \endverbatim
171*>
172*> \param[out] INFO
173*> \verbatim
174*> INFO is INTEGER
175*> = 0: successful exit
176*> < 0: if INFO = -i, the i-th argument had an illegal value
177*> > 0: DPOTRF or DSYEVD returned an error code:
178*> <= N: if INFO = i and JOBZ = 'N', then the algorithm
179*> failed to converge; i off-diagonal elements of an
180*> intermediate tridiagonal form did not converge to
181*> zero;
182*> if INFO = i and JOBZ = 'V', then the algorithm
183*> failed to compute an eigenvalue while working on
184*> the submatrix lying in rows and columns INFO/(N+1)
185*> through mod(INFO,N+1);
186*> > N: if INFO = N + i, for 1 <= i <= N, then the leading
187*> principal minor of order i of B is not positive.
188*> The factorization of B could not be completed and
189*> no eigenvalues or eigenvectors were computed.
190*> \endverbatim
191*
192* Authors:
193* ========
194*
195*> \author Univ. of Tennessee
196*> \author Univ. of California Berkeley
197*> \author Univ. of Colorado Denver
198*> \author NAG Ltd.
199*
200*> \ingroup hegvd
201*
202*> \par Further Details:
203* =====================
204*>
205*> \verbatim
206*>
207*> Modified so that no backsubstitution is performed if DSYEVD fails to
208*> converge (NEIG in old code could be greater than N causing out of
209*> bounds reference to A - reported by Ralf Meyer). Also corrected the
210*> description of INFO and the test on ITYPE. Sven, 16 Feb 05.
211*> \endverbatim
212*
213*> \par Contributors:
214* ==================
215*>
216*> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
217*>
218* =====================================================================
219 SUBROUTINE dsygvd( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
220 $ LWORK, IWORK, LIWORK, INFO )
221*
222* -- LAPACK driver routine --
223* -- LAPACK is a software package provided by Univ. of Tennessee, --
224* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225*
226* .. Scalar Arguments ..
227 CHARACTER JOBZ, UPLO
228 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
229* ..
230* .. Array Arguments ..
231 INTEGER IWORK( * )
232 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
233* ..
234*
235* =====================================================================
236*
237* .. Parameters ..
238 DOUBLE PRECISION ONE
239 parameter( one = 1.0d+0 )
240* ..
241* .. Local Scalars ..
242 LOGICAL LQUERY, UPPER, WANTZ
243 CHARACTER TRANS
244 INTEGER LIOPT, LIWMIN, LOPT, LWMIN
245* ..
246* .. External Functions ..
247 LOGICAL LSAME
248 EXTERNAL lsame
249* ..
250* .. External Subroutines ..
251 EXTERNAL dpotrf, dsyevd, dsygst, dtrmm, dtrsm, xerbla
252* ..
253* .. Intrinsic Functions ..
254 INTRINSIC dble, max
255* ..
256* .. Executable Statements ..
257*
258* Test the input parameters.
259*
260 wantz = lsame( jobz, 'V' )
261 upper = lsame( uplo, 'U' )
262 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
263*
264 info = 0
265 IF( n.LE.1 ) THEN
266 liwmin = 1
267 lwmin = 1
268 ELSE IF( wantz ) THEN
269 liwmin = 3 + 5*n
270 lwmin = 1 + 6*n + 2*n**2
271 ELSE
272 liwmin = 1
273 lwmin = 2*n + 1
274 END IF
275 lopt = lwmin
276 liopt = liwmin
277 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
278 info = -1
279 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
280 info = -2
281 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
282 info = -3
283 ELSE IF( n.LT.0 ) THEN
284 info = -4
285 ELSE IF( lda.LT.max( 1, n ) ) THEN
286 info = -6
287 ELSE IF( ldb.LT.max( 1, n ) ) THEN
288 info = -8
289 END IF
290*
291 IF( info.EQ.0 ) THEN
292 work( 1 ) = lopt
293 iwork( 1 ) = liopt
294*
295 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
296 info = -11
297 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
298 info = -13
299 END IF
300 END IF
301*
302 IF( info.NE.0 ) THEN
303 CALL xerbla( 'DSYGVD', -info )
304 RETURN
305 ELSE IF( lquery ) THEN
306 RETURN
307 END IF
308*
309* Quick return if possible
310*
311 IF( n.EQ.0 )
312 $ RETURN
313*
314* Form a Cholesky factorization of B.
315*
316 CALL dpotrf( uplo, n, b, ldb, info )
317 IF( info.NE.0 ) THEN
318 info = n + info
319 RETURN
320 END IF
321*
322* Transform problem to standard eigenvalue problem and solve.
323*
324 CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
325 CALL dsyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork,
326 $ info )
327 lopt = int( max( dble( lopt ), dble( work( 1 ) ) ) )
328 liopt = int( max( dble( liopt ), dble( iwork( 1 ) ) ) )
329*
330 IF( wantz .AND. info.EQ.0 ) THEN
331*
332* Backtransform eigenvectors to the original problem.
333*
334 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
335*
336* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
337* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
338*
339 IF( upper ) THEN
340 trans = 'N'
341 ELSE
342 trans = 'T'
343 END IF
344*
345 CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, n, one,
346 $ b, ldb, a, lda )
347*
348 ELSE IF( itype.EQ.3 ) THEN
349*
350* For B*A*x=(lambda)*x;
351* backtransform eigenvectors: x = L*y or U**T*y
352*
353 IF( upper ) THEN
354 trans = 'T'
355 ELSE
356 trans = 'N'
357 END IF
358*
359 CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, n, one,
360 $ b, ldb, a, lda )
361 END IF
362 END IF
363*
364 work( 1 ) = lopt
365 iwork( 1 ) = liopt
366*
367 RETURN
368*
369* End of DSYGVD
370*
371 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
DSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition dsyevd.f:178
subroutine dsygst(itype, uplo, n, a, lda, b, ldb, info)
DSYGST
Definition dsygst.f:127
subroutine dsygvd(itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, iwork, liwork, info)
DSYGVD
Definition dsygvd.f:221
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181