LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sspevd.f
Go to the documentation of this file.
1*> \brief <b> SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SSPEVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sspevd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sspevd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sspevd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
22* IWORK, LIWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBZ, UPLO
26* INTEGER INFO, LDZ, LIWORK, LWORK, N
27* ..
28* .. Array Arguments ..
29* INTEGER IWORK( * )
30* REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SSPEVD computes all the eigenvalues and, optionally, eigenvectors
40*> of a real symmetric matrix A in packed storage. If eigenvectors are
41*> desired, it uses a divide and conquer algorithm.
42*>
43*> The divide and conquer algorithm makes very mild assumptions about
44*> floating point arithmetic. It will work on machines with a guard
45*> digit in add/subtract, or on those binary machines without guard
46*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
47*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
48*> without guard digits, but we know of none.
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] JOBZ
55*> \verbatim
56*> JOBZ is CHARACTER*1
57*> = 'N': Compute eigenvalues only;
58*> = 'V': Compute eigenvalues and eigenvectors.
59*> \endverbatim
60*>
61*> \param[in] UPLO
62*> \verbatim
63*> UPLO is CHARACTER*1
64*> = 'U': Upper triangle of A is stored;
65*> = 'L': Lower triangle of A is stored.
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*> N is INTEGER
71*> The order of the matrix A. N >= 0.
72*> \endverbatim
73*>
74*> \param[in,out] AP
75*> \verbatim
76*> AP is REAL array, dimension (N*(N+1)/2)
77*> On entry, the upper or lower triangle of the symmetric matrix
78*> A, packed columnwise in a linear array. The j-th column of A
79*> is stored in the array AP as follows:
80*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
81*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
82*>
83*> On exit, AP is overwritten by values generated during the
84*> reduction to tridiagonal form. If UPLO = 'U', the diagonal
85*> and first superdiagonal of the tridiagonal matrix T overwrite
86*> the corresponding elements of A, and if UPLO = 'L', the
87*> diagonal and first subdiagonal of T overwrite the
88*> corresponding elements of A.
89*> \endverbatim
90*>
91*> \param[out] W
92*> \verbatim
93*> W is REAL array, dimension (N)
94*> If INFO = 0, the eigenvalues in ascending order.
95*> \endverbatim
96*>
97*> \param[out] Z
98*> \verbatim
99*> Z is REAL array, dimension (LDZ, N)
100*> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
101*> eigenvectors of the matrix A, with the i-th column of Z
102*> holding the eigenvector associated with W(i).
103*> If JOBZ = 'N', then Z is not referenced.
104*> \endverbatim
105*>
106*> \param[in] LDZ
107*> \verbatim
108*> LDZ is INTEGER
109*> The leading dimension of the array Z. LDZ >= 1, and if
110*> JOBZ = 'V', LDZ >= max(1,N).
111*> \endverbatim
112*>
113*> \param[out] WORK
114*> \verbatim
115*> WORK is REAL array, dimension (MAX(1,LWORK))
116*> On exit, if INFO = 0, WORK(1) returns the required LWORK.
117*> \endverbatim
118*>
119*> \param[in] LWORK
120*> \verbatim
121*> LWORK is INTEGER
122*> The dimension of the array WORK.
123*> If N <= 1, LWORK must be at least 1.
124*> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
125*> If JOBZ = 'V' and N > 1, LWORK must be at least
126*> 1 + 6*N + N**2.
127*>
128*> If LWORK = -1, then a workspace query is assumed; the routine
129*> only calculates the required sizes of the WORK and IWORK
130*> arrays, returns these values as the first entries of the WORK
131*> and IWORK arrays, and no error message related to LWORK or
132*> LIWORK is issued by XERBLA.
133*> \endverbatim
134*>
135*> \param[out] IWORK
136*> \verbatim
137*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
138*> On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
139*> \endverbatim
140*>
141*> \param[in] LIWORK
142*> \verbatim
143*> LIWORK is INTEGER
144*> The dimension of the array IWORK.
145*> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
146*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
147*>
148*> If LIWORK = -1, then a workspace query is assumed; the
149*> routine only calculates the required sizes of the WORK and
150*> IWORK arrays, returns these values as the first entries of
151*> the WORK and IWORK arrays, and no error message related to
152*> LWORK or LIWORK is issued by XERBLA.
153*> \endverbatim
154*>
155*> \param[out] INFO
156*> \verbatim
157*> INFO is INTEGER
158*> = 0: successful exit
159*> < 0: if INFO = -i, the i-th argument had an illegal value.
160*> > 0: if INFO = i, the algorithm failed to converge; i
161*> off-diagonal elements of an intermediate tridiagonal
162*> form did not converge to zero.
163*> \endverbatim
164*
165* Authors:
166* ========
167*
168*> \author Univ. of Tennessee
169*> \author Univ. of California Berkeley
170*> \author Univ. of Colorado Denver
171*> \author NAG Ltd.
172*
173*> \ingroup realOTHEReigen
174*
175* =====================================================================
176 SUBROUTINE sspevd( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
177 $ IWORK, LIWORK, INFO )
178*
179* -- LAPACK driver routine --
180* -- LAPACK is a software package provided by Univ. of Tennessee, --
181* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182*
183* .. Scalar Arguments ..
184 CHARACTER JOBZ, UPLO
185 INTEGER INFO, LDZ, LIWORK, LWORK, N
186* ..
187* .. Array Arguments ..
188 INTEGER IWORK( * )
189 REAL AP( * ), W( * ), WORK( * ), Z( LDZ, * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 REAL ZERO, ONE
196 parameter( zero = 0.0e+0, one = 1.0e+0 )
197* ..
198* .. Local Scalars ..
199 LOGICAL LQUERY, WANTZ
200 INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
201 $ llwork, lwmin
202 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
203 $ smlnum
204* ..
205* .. External Functions ..
206 LOGICAL LSAME
207 REAL SLAMCH, SLANSP
208 EXTERNAL lsame, slamch, slansp
209* ..
210* .. External Subroutines ..
211 EXTERNAL sopmtr, sscal, ssptrd, sstedc, ssterf, xerbla
212* ..
213* .. Intrinsic Functions ..
214 INTRINSIC sqrt
215* ..
216* .. Executable Statements ..
217*
218* Test the input parameters.
219*
220 wantz = lsame( jobz, 'V' )
221 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
222*
223 info = 0
224 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
225 info = -1
226 ELSE IF( .NOT.( lsame( uplo, 'U' ) .OR. lsame( uplo, 'L' ) ) )
227 $ THEN
228 info = -2
229 ELSE IF( n.LT.0 ) THEN
230 info = -3
231 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
232 info = -7
233 END IF
234*
235 IF( info.EQ.0 ) THEN
236 IF( n.LE.1 ) THEN
237 liwmin = 1
238 lwmin = 1
239 ELSE
240 IF( wantz ) THEN
241 liwmin = 3 + 5*n
242 lwmin = 1 + 6*n + n**2
243 ELSE
244 liwmin = 1
245 lwmin = 2*n
246 END IF
247 END IF
248 iwork( 1 ) = liwmin
249 work( 1 ) = lwmin
250*
251 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
252 info = -9
253 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
254 info = -11
255 END IF
256 END IF
257*
258 IF( info.NE.0 ) THEN
259 CALL xerbla( 'SSPEVD', -info )
260 RETURN
261 ELSE IF( lquery ) THEN
262 RETURN
263 END IF
264*
265* Quick return if possible
266*
267 IF( n.EQ.0 )
268 $ RETURN
269*
270 IF( n.EQ.1 ) THEN
271 w( 1 ) = ap( 1 )
272 IF( wantz )
273 $ z( 1, 1 ) = one
274 RETURN
275 END IF
276*
277* Get machine constants.
278*
279 safmin = slamch( 'Safe minimum' )
280 eps = slamch( 'Precision' )
281 smlnum = safmin / eps
282 bignum = one / smlnum
283 rmin = sqrt( smlnum )
284 rmax = sqrt( bignum )
285*
286* Scale matrix to allowable range, if necessary.
287*
288 anrm = slansp( 'M', uplo, n, ap, work )
289 iscale = 0
290 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
291 iscale = 1
292 sigma = rmin / anrm
293 ELSE IF( anrm.GT.rmax ) THEN
294 iscale = 1
295 sigma = rmax / anrm
296 END IF
297 IF( iscale.EQ.1 ) THEN
298 CALL sscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
299 END IF
300*
301* Call SSPTRD to reduce symmetric packed matrix to tridiagonal form.
302*
303 inde = 1
304 indtau = inde + n
305 CALL ssptrd( uplo, n, ap, w, work( inde ), work( indtau ), iinfo )
306*
307* For eigenvalues only, call SSTERF. For eigenvectors, first call
308* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
309* tridiagonal matrix, then call SOPMTR to multiply it by the
310* Householder transformations represented in AP.
311*
312 IF( .NOT.wantz ) THEN
313 CALL ssterf( n, w, work( inde ), info )
314 ELSE
315 indwrk = indtau + n
316 llwork = lwork - indwrk + 1
317 CALL sstedc( 'I', n, w, work( inde ), z, ldz, work( indwrk ),
318 $ llwork, iwork, liwork, info )
319 CALL sopmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
320 $ work( indwrk ), iinfo )
321 END IF
322*
323* If matrix was scaled, then rescale eigenvalues appropriately.
324*
325 IF( iscale.EQ.1 )
326 $ CALL sscal( n, one / sigma, w, 1 )
327*
328 work( 1 ) = lwmin
329 iwork( 1 ) = liwmin
330 RETURN
331*
332* End of SSPEVD
333*
334 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:188
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine ssptrd(UPLO, N, AP, D, E, TAU, INFO)
SSPTRD
Definition: ssptrd.f:150
subroutine sopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
SOPMTR
Definition: sopmtr.f:150
subroutine sspevd(JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: sspevd.f:178
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79