LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssygvd.f
Go to the documentation of this file.
1*> \brief \b SSYGVD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SSYGVD + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssygvd.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssygvd.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssygvd.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
20* LWORK, IWORK, LIWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER JOBZ, UPLO
24* INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
25* ..
26* .. Array Arguments ..
27* INTEGER IWORK( * )
28* REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SSYGVD computes all the eigenvalues, and optionally, the eigenvectors
38*> of a real generalized symmetric-definite eigenproblem, of the form
39*> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and
40*> B are assumed to be symmetric and B is also positive definite.
41*> If eigenvectors are desired, it uses a divide and conquer algorithm.
42*>
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 REAL 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 REAL array, dimension (LDB, N)
105*> On entry, the symmetric matrix B. If UPLO = 'U', the
106*> leading N-by-N upper triangular part of B contains the
107*> upper triangular part of the matrix B. If UPLO = 'L',
108*> the leading N-by-N lower triangular part of B contains
109*> 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 REAL array, dimension (N)
125*> If INFO = 0, the eigenvalues in ascending order.
126*> \endverbatim
127*>
128*> \param[out] WORK
129*> \verbatim
130*> WORK is REAL 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 dimension of the array WORK.
138*> If N <= 1, LWORK >= 1.
139*> If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
140*> If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
141*>
142*> If LWORK = -1, then a workspace query is assumed; the routine
143*> only calculates the optimal sizes of the WORK and IWORK
144*> arrays, returns these values as the first entries of the WORK
145*> and IWORK arrays, and no error message related to LWORK or
146*> LIWORK is issued by XERBLA.
147*> \endverbatim
148*>
149*> \param[out] IWORK
150*> \verbatim
151*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
152*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
153*> \endverbatim
154*>
155*> \param[in] LIWORK
156*> \verbatim
157*> LIWORK is INTEGER
158*> The dimension of the array IWORK.
159*> If N <= 1, LIWORK >= 1.
160*> If JOBZ = 'N' and N > 1, LIWORK >= 1.
161*> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
162*>
163*> If LIWORK = -1, then a workspace query is assumed; the
164*> routine only calculates the optimal sizes of the WORK and
165*> IWORK arrays, returns these values as the first entries of
166*> the WORK and IWORK arrays, and no error message related to
167*> LWORK or LIWORK is issued by XERBLA.
168*> \endverbatim
169*>
170*> \param[out] INFO
171*> \verbatim
172*> INFO is INTEGER
173*> = 0: successful exit
174*> < 0: if INFO = -i, the i-th argument had an illegal value
175*> > 0: SPOTRF or SSYEVD returned an error code:
176*> <= N: if INFO = i and JOBZ = 'N', then the algorithm
177*> failed to converge; i off-diagonal elements of an
178*> intermediate tridiagonal form did not converge to
179*> zero;
180*> if INFO = i and JOBZ = 'V', then the algorithm
181*> failed to compute an eigenvalue while working on
182*> the submatrix lying in rows and columns INFO/(N+1)
183*> through mod(INFO,N+1);
184*> > N: if INFO = N + i, for 1 <= i <= N, then the leading
185*> principal minor of order i of B is not positive.
186*> The factorization of B could not be completed and
187*> no eigenvalues or eigenvectors were computed.
188*> \endverbatim
189*
190* Authors:
191* ========
192*
193*> \author Univ. of Tennessee
194*> \author Univ. of California Berkeley
195*> \author Univ. of Colorado Denver
196*> \author NAG Ltd.
197*
198*> \ingroup hegvd
199*
200*> \par Further Details:
201* =====================
202*>
203*> \verbatim
204*>
205*> Modified so that no backsubstitution is performed if SSYEVD fails to
206*> converge (NEIG in old code could be greater than N causing out of
207*> bounds reference to A - reported by Ralf Meyer). Also corrected the
208*> description of INFO and the test on ITYPE. Sven, 16 Feb 05.
209*> \endverbatim
210*
211*> \par Contributors:
212* ==================
213*>
214*> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
215*>
216* =====================================================================
217 SUBROUTINE ssygvd( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W,
218 $ WORK,
219 $ LWORK, IWORK, LIWORK, INFO )
220*
221* -- LAPACK driver routine --
222* -- LAPACK is a software package provided by Univ. of Tennessee, --
223* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224*
225* .. Scalar Arguments ..
226 CHARACTER JOBZ, UPLO
227 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
228* ..
229* .. Array Arguments ..
230 INTEGER IWORK( * )
231 REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
232* ..
233*
234* =====================================================================
235*
236* .. Parameters ..
237 REAL ONE
238 PARAMETER ( ONE = 1.0e+0 )
239* ..
240* .. Local Scalars ..
241 LOGICAL LQUERY, UPPER, WANTZ
242 CHARACTER TRANS
243 INTEGER LIOPT, LIWMIN, LOPT, LWMIN
244* ..
245* .. External Functions ..
246 LOGICAL LSAME
247 REAL SROUNDUP_LWORK
248 EXTERNAL lsame, sroundup_lwork
249* ..
250* .. External Subroutines ..
251 EXTERNAL spotrf, ssyevd, ssygst, strmm, strsm,
252 $ xerbla
253* ..
254* .. Intrinsic Functions ..
255 INTRINSIC max, real
256* ..
257* .. Executable Statements ..
258*
259* Test the input parameters.
260*
261 wantz = lsame( jobz, 'V' )
262 upper = lsame( uplo, 'U' )
263 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
264*
265 info = 0
266 IF( n.LE.1 ) THEN
267 liwmin = 1
268 lwmin = 1
269 ELSE IF( wantz ) THEN
270 liwmin = 3 + 5*n
271 lwmin = 1 + 6*n + 2*n**2
272 ELSE
273 liwmin = 1
274 lwmin = 2*n + 1
275 END IF
276 lopt = lwmin
277 liopt = liwmin
278 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
279 info = -1
280 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
281 info = -2
282 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
283 info = -3
284 ELSE IF( n.LT.0 ) THEN
285 info = -4
286 ELSE IF( lda.LT.max( 1, n ) ) THEN
287 info = -6
288 ELSE IF( ldb.LT.max( 1, n ) ) THEN
289 info = -8
290 END IF
291*
292 IF( info.EQ.0 ) THEN
293 work( 1 ) = sroundup_lwork(lopt)
294 iwork( 1 ) = liopt
295*
296 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
297 info = -11
298 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
299 info = -13
300 END IF
301 END IF
302*
303 IF( info.NE.0 ) THEN
304 CALL xerbla( 'SSYGVD', -info )
305 RETURN
306 ELSE IF( lquery ) THEN
307 RETURN
308 END IF
309*
310* Quick return if possible
311*
312 IF( n.EQ.0 )
313 $ RETURN
314*
315* Form a Cholesky factorization of B.
316*
317 CALL spotrf( uplo, n, b, ldb, info )
318 IF( info.NE.0 ) THEN
319 info = n + info
320 RETURN
321 END IF
322*
323* Transform problem to standard eigenvalue problem and solve.
324*
325 CALL ssygst( itype, uplo, n, a, lda, b, ldb, info )
326 CALL ssyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork,
327 $ liwork,
328 $ info )
329 lopt = int( max( real( lopt ), real( work( 1 ) ) ) )
330 liopt = int( max( real( liopt ), real( iwork( 1 ) ) ) )
331*
332 IF( wantz .AND. info.EQ.0 ) THEN
333*
334* Backtransform eigenvectors to the original problem.
335*
336 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
337*
338* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
339* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
340*
341 IF( upper ) THEN
342 trans = 'N'
343 ELSE
344 trans = 'T'
345 END IF
346*
347 CALL strsm( 'Left', uplo, trans, 'Non-unit', n, n, one,
348 $ b, ldb, a, lda )
349*
350 ELSE IF( itype.EQ.3 ) THEN
351*
352* For B*A*x=(lambda)*x;
353* backtransform eigenvectors: x = L*y or U**T*y
354*
355 IF( upper ) THEN
356 trans = 'T'
357 ELSE
358 trans = 'N'
359 END IF
360*
361 CALL strmm( 'Left', uplo, trans, 'Non-unit', n, n, one,
362 $ b, ldb, a, lda )
363 END IF
364 END IF
365*
366 work( 1 ) = sroundup_lwork(lopt)
367 iwork( 1 ) = liopt
368*
369 RETURN
370*
371* End of SSYGVD
372*
373 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ssyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
SSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition ssyevd.f:174
subroutine ssygst(itype, uplo, n, a, lda, b, ldb, info)
SSYGST
Definition ssygst.f:125
subroutine ssygvd(itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, iwork, liwork, info)
SSYGVD
Definition ssygvd.f:220
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:105
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181