LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dsygv.f
Go to the documentation of this file.
1 *> \brief \b DSYGST
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSYGV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
22 * LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DSYGV computes all the eigenvalues, and optionally, the eigenvectors
39 *> of a real generalized symmetric-definite eigenproblem, of the form
40 *> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
41 *> Here A and B are assumed to be symmetric and B is also
42 *> positive definite.
43 *> \endverbatim
44 *
45 * Arguments:
46 * ==========
47 *
48 *> \param[in] ITYPE
49 *> \verbatim
50 *> ITYPE is INTEGER
51 *> Specifies the problem type to be solved:
52 *> = 1: A*x = (lambda)*B*x
53 *> = 2: A*B*x = (lambda)*x
54 *> = 3: B*A*x = (lambda)*x
55 *> \endverbatim
56 *>
57 *> \param[in] JOBZ
58 *> \verbatim
59 *> JOBZ is CHARACTER*1
60 *> = 'N': Compute eigenvalues only;
61 *> = 'V': Compute eigenvalues and eigenvectors.
62 *> \endverbatim
63 *>
64 *> \param[in] UPLO
65 *> \verbatim
66 *> UPLO is CHARACTER*1
67 *> = 'U': Upper triangles of A and B are stored;
68 *> = 'L': Lower triangles of A and B are stored.
69 *> \endverbatim
70 *>
71 *> \param[in] N
72 *> \verbatim
73 *> N is INTEGER
74 *> The order of the matrices A and B. N >= 0.
75 *> \endverbatim
76 *>
77 *> \param[in,out] A
78 *> \verbatim
79 *> A is DOUBLE PRECISION array, dimension (LDA, N)
80 *> On entry, the symmetric matrix A. If UPLO = 'U', the
81 *> leading N-by-N upper triangular part of A contains the
82 *> upper triangular part of the matrix A. If UPLO = 'L',
83 *> the leading N-by-N lower triangular part of A contains
84 *> the lower triangular part of the matrix A.
85 *>
86 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
87 *> matrix Z of eigenvectors. The eigenvectors are normalized
88 *> as follows:
89 *> if ITYPE = 1 or 2, Z**T*B*Z = I;
90 *> if ITYPE = 3, Z**T*inv(B)*Z = I.
91 *> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
92 *> or the lower triangle (if UPLO='L') of A, including the
93 *> diagonal, is destroyed.
94 *> \endverbatim
95 *>
96 *> \param[in] LDA
97 *> \verbatim
98 *> LDA is INTEGER
99 *> The leading dimension of the array A. LDA >= max(1,N).
100 *> \endverbatim
101 *>
102 *> \param[in,out] B
103 *> \verbatim
104 *> B is DOUBLE PRECISION array, dimension (LDB, N)
105 *> On entry, the symmetric positive definite matrix B.
106 *> If UPLO = 'U', the leading N-by-N upper triangular part of B
107 *> contains the upper triangular part of the matrix B.
108 *> If UPLO = 'L', the leading N-by-N lower triangular part of B
109 *> contains the lower triangular part of the matrix B.
110 *>
111 *> On exit, if INFO <= N, the part of B containing the matrix is
112 *> overwritten by the triangular factor U or L from the Cholesky
113 *> factorization B = U**T*U or B = L*L**T.
114 *> \endverbatim
115 *>
116 *> \param[in] LDB
117 *> \verbatim
118 *> LDB is INTEGER
119 *> The leading dimension of the array B. LDB >= max(1,N).
120 *> \endverbatim
121 *>
122 *> \param[out] W
123 *> \verbatim
124 *> W is DOUBLE PRECISION array, dimension (N)
125 *> If INFO = 0, the eigenvalues in ascending order.
126 *> \endverbatim
127 *>
128 *> \param[out] WORK
129 *> \verbatim
130 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
131 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
132 *> \endverbatim
133 *>
134 *> \param[in] LWORK
135 *> \verbatim
136 *> LWORK is INTEGER
137 *> The length of the array WORK. LWORK >= max(1,3*N-1).
138 *> For optimal efficiency, LWORK >= (NB+2)*N,
139 *> where NB is the blocksize for DSYTRD returned by ILAENV.
140 *>
141 *> If LWORK = -1, then a workspace query is assumed; the routine
142 *> only calculates the optimal size of the WORK array, returns
143 *> this value as the first entry of the WORK array, and no error
144 *> message related to LWORK is issued by XERBLA.
145 *> \endverbatim
146 *>
147 *> \param[out] INFO
148 *> \verbatim
149 *> INFO is INTEGER
150 *> = 0: successful exit
151 *> < 0: if INFO = -i, the i-th argument had an illegal value
152 *> > 0: DPOTRF or DSYEV returned an error code:
153 *> <= N: if INFO = i, DSYEV failed to converge;
154 *> i off-diagonal elements of an intermediate
155 *> tridiagonal form did not converge to zero;
156 *> > N: if INFO = N + i, for 1 <= i <= N, then the leading
157 *> minor of order i of B is not positive definite.
158 *> The factorization of B could not be completed and
159 *> no eigenvalues or eigenvectors were computed.
160 *> \endverbatim
161 *
162 * Authors:
163 * ========
164 *
165 *> \author Univ. of Tennessee
166 *> \author Univ. of California Berkeley
167 *> \author Univ. of Colorado Denver
168 *> \author NAG Ltd.
169 *
170 *> \date November 2011
171 *
172 *> \ingroup doubleSYeigen
173 *
174 * =====================================================================
175  SUBROUTINE dsygv( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
176  $ lwork, info )
177 *
178 * -- LAPACK driver routine (version 3.4.0) --
179 * -- LAPACK is a software package provided by Univ. of Tennessee, --
180 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181 * November 2011
182 *
183 * .. Scalar Arguments ..
184  CHARACTER jobz, uplo
185  INTEGER info, itype, lda, ldb, lwork, n
186 * ..
187 * .. Array Arguments ..
188  DOUBLE PRECISION a( lda, * ), b( ldb, * ), w( * ), work( * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  DOUBLE PRECISION one
195  parameter( one = 1.0d+0 )
196 * ..
197 * .. Local Scalars ..
198  LOGICAL lquery, upper, wantz
199  CHARACTER trans
200  INTEGER lwkmin, lwkopt, nb, neig
201 * ..
202 * .. External Functions ..
203  LOGICAL lsame
204  INTEGER ilaenv
205  EXTERNAL lsame, ilaenv
206 * ..
207 * .. External Subroutines ..
208  EXTERNAL dpotrf, dsyev, dsygst, dtrmm, dtrsm, xerbla
209 * ..
210 * .. Intrinsic Functions ..
211  INTRINSIC max
212 * ..
213 * .. Executable Statements ..
214 *
215 * Test the input parameters.
216 *
217  wantz = lsame( jobz, 'V' )
218  upper = lsame( uplo, 'U' )
219  lquery = ( lwork.EQ.-1 )
220 *
221  info = 0
222  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
223  info = -1
224  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
225  info = -2
226  ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
227  info = -3
228  ELSE IF( n.LT.0 ) THEN
229  info = -4
230  ELSE IF( lda.LT.max( 1, n ) ) THEN
231  info = -6
232  ELSE IF( ldb.LT.max( 1, n ) ) THEN
233  info = -8
234  END IF
235 *
236  IF( info.EQ.0 ) THEN
237  lwkmin = max( 1, 3*n - 1 )
238  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
239  lwkopt = max( lwkmin, ( nb + 2 )*n )
240  work( 1 ) = lwkopt
241 *
242  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
243  info = -11
244  END IF
245  END IF
246 *
247  IF( info.NE.0 ) THEN
248  CALL xerbla( 'DSYGV ', -info )
249  return
250  ELSE IF( lquery ) THEN
251  return
252  END IF
253 *
254 * Quick return if possible
255 *
256  IF( n.EQ.0 )
257  $ return
258 *
259 * Form a Cholesky factorization of B.
260 *
261  CALL dpotrf( uplo, n, b, ldb, info )
262  IF( info.NE.0 ) THEN
263  info = n + info
264  return
265  END IF
266 *
267 * Transform problem to standard eigenvalue problem and solve.
268 *
269  CALL dsygst( itype, uplo, n, a, lda, b, ldb, info )
270  CALL dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )
271 *
272  IF( wantz ) THEN
273 *
274 * Backtransform eigenvectors to the original problem.
275 *
276  neig = n
277  IF( info.GT.0 )
278  $ neig = info - 1
279  IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
280 *
281 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
282 * backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
283 *
284  IF( upper ) THEN
285  trans = 'N'
286  ELSE
287  trans = 'T'
288  END IF
289 *
290  CALL dtrsm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
291  $ b, ldb, a, lda )
292 *
293  ELSE IF( itype.EQ.3 ) THEN
294 *
295 * For B*A*x=(lambda)*x;
296 * backtransform eigenvectors: x = L*y or U**T*y
297 *
298  IF( upper ) THEN
299  trans = 'T'
300  ELSE
301  trans = 'N'
302  END IF
303 *
304  CALL dtrmm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
305  $ b, ldb, a, lda )
306  END IF
307  END IF
308 *
309  work( 1 ) = lwkopt
310  return
311 *
312 * End of DSYGV
313 *
314  END