LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sstevx.f
Go to the documentation of this file.
1*> \brief <b> SSTEVX 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 SSTEVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstevx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstevx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstevx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
22* M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBZ, RANGE
26* INTEGER IL, INFO, IU, LDZ, M, N
27* REAL ABSTOL, VL, VU
28* ..
29* .. Array Arguments ..
30* INTEGER IFAIL( * ), IWORK( * )
31* REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> SSTEVX computes selected eigenvalues and, optionally, eigenvectors
41*> of a real symmetric tridiagonal matrix A. Eigenvalues and
42*> eigenvectors can be selected by specifying either a range of values
43*> or a range of indices for the desired eigenvalues.
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] JOBZ
50*> \verbatim
51*> JOBZ is CHARACTER*1
52*> = 'N': Compute eigenvalues only;
53*> = 'V': Compute eigenvalues and eigenvectors.
54*> \endverbatim
55*>
56*> \param[in] RANGE
57*> \verbatim
58*> RANGE is CHARACTER*1
59*> = 'A': all eigenvalues will be found.
60*> = 'V': all eigenvalues in the half-open interval (VL,VU]
61*> will be found.
62*> = 'I': the IL-th through IU-th eigenvalues will be found.
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*> N is INTEGER
68*> The order of the matrix. N >= 0.
69*> \endverbatim
70*>
71*> \param[in,out] D
72*> \verbatim
73*> D is REAL array, dimension (N)
74*> On entry, the n diagonal elements of the tridiagonal matrix
75*> A.
76*> On exit, D may be multiplied by a constant factor chosen
77*> to avoid over/underflow in computing the eigenvalues.
78*> \endverbatim
79*>
80*> \param[in,out] E
81*> \verbatim
82*> E is REAL array, dimension (max(1,N-1))
83*> On entry, the (n-1) subdiagonal elements of the tridiagonal
84*> matrix A in elements 1 to N-1 of E.
85*> On exit, E may be multiplied by a constant factor chosen
86*> to avoid over/underflow in computing the eigenvalues.
87*> \endverbatim
88*>
89*> \param[in] VL
90*> \verbatim
91*> VL is REAL
92*> If RANGE='V', the lower bound of the interval to
93*> be searched for eigenvalues. VL < VU.
94*> Not referenced if RANGE = 'A' or 'I'.
95*> \endverbatim
96*>
97*> \param[in] VU
98*> \verbatim
99*> VU is REAL
100*> If RANGE='V', the upper bound of the interval to
101*> be searched for eigenvalues. VL < VU.
102*> Not referenced if RANGE = 'A' or 'I'.
103*> \endverbatim
104*>
105*> \param[in] IL
106*> \verbatim
107*> IL is INTEGER
108*> If RANGE='I', the index of the
109*> smallest eigenvalue to be returned.
110*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
111*> Not referenced if RANGE = 'A' or 'V'.
112*> \endverbatim
113*>
114*> \param[in] IU
115*> \verbatim
116*> IU is INTEGER
117*> If RANGE='I', the index of the
118*> largest eigenvalue to be returned.
119*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
120*> Not referenced if RANGE = 'A' or 'V'.
121*> \endverbatim
122*>
123*> \param[in] ABSTOL
124*> \verbatim
125*> ABSTOL is REAL
126*> The absolute error tolerance for the eigenvalues.
127*> An approximate eigenvalue is accepted as converged
128*> when it is determined to lie in an interval [a,b]
129*> of width less than or equal to
130*>
131*> ABSTOL + EPS * max( |a|,|b| ) ,
132*>
133*> where EPS is the machine precision. If ABSTOL is less
134*> than or equal to zero, then EPS*|T| will be used in
135*> its place, where |T| is the 1-norm of the tridiagonal
136*> matrix.
137*>
138*> Eigenvalues will be computed most accurately when ABSTOL is
139*> set to twice the underflow threshold 2*SLAMCH('S'), not zero.
140*> If this routine returns with INFO>0, indicating that some
141*> eigenvectors did not converge, try setting ABSTOL to
142*> 2*SLAMCH('S').
143*>
144*> See "Computing Small Singular Values of Bidiagonal Matrices
145*> with Guaranteed High Relative Accuracy," by Demmel and
146*> Kahan, LAPACK Working Note #3.
147*> \endverbatim
148*>
149*> \param[out] M
150*> \verbatim
151*> M is INTEGER
152*> The total number of eigenvalues found. 0 <= M <= N.
153*> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
154*> \endverbatim
155*>
156*> \param[out] W
157*> \verbatim
158*> W is REAL array, dimension (N)
159*> The first M elements contain the selected eigenvalues in
160*> ascending order.
161*> \endverbatim
162*>
163*> \param[out] Z
164*> \verbatim
165*> Z is REAL array, dimension (LDZ, max(1,M) )
166*> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
167*> contain the orthonormal eigenvectors of the matrix A
168*> corresponding to the selected eigenvalues, with the i-th
169*> column of Z holding the eigenvector associated with W(i).
170*> If an eigenvector fails to converge (INFO > 0), then that
171*> column of Z contains the latest approximation to the
172*> eigenvector, and the index of the eigenvector is returned
173*> in IFAIL. If JOBZ = 'N', then Z is not referenced.
174*> Note: the user must ensure that at least max(1,M) columns are
175*> supplied in the array Z; if RANGE = 'V', the exact value of M
176*> is not known in advance and an upper bound must be used.
177*> \endverbatim
178*>
179*> \param[in] LDZ
180*> \verbatim
181*> LDZ is INTEGER
182*> The leading dimension of the array Z. LDZ >= 1, and if
183*> JOBZ = 'V', LDZ >= max(1,N).
184*> \endverbatim
185*>
186*> \param[out] WORK
187*> \verbatim
188*> WORK is REAL array, dimension (5*N)
189*> \endverbatim
190*>
191*> \param[out] IWORK
192*> \verbatim
193*> IWORK is INTEGER array, dimension (5*N)
194*> \endverbatim
195*>
196*> \param[out] IFAIL
197*> \verbatim
198*> IFAIL is INTEGER array, dimension (N)
199*> If JOBZ = 'V', then if INFO = 0, the first M elements of
200*> IFAIL are zero. If INFO > 0, then IFAIL contains the
201*> indices of the eigenvectors that failed to converge.
202*> If JOBZ = 'N', then IFAIL is not referenced.
203*> \endverbatim
204*>
205*> \param[out] INFO
206*> \verbatim
207*> INFO is INTEGER
208*> = 0: successful exit
209*> < 0: if INFO = -i, the i-th argument had an illegal value
210*> > 0: if INFO = i, then i eigenvectors failed to converge.
211*> Their indices are stored in array IFAIL.
212*> \endverbatim
213*
214* Authors:
215* ========
216*
217*> \author Univ. of Tennessee
218*> \author Univ. of California Berkeley
219*> \author Univ. of Colorado Denver
220*> \author NAG Ltd.
221*
222*> \ingroup realOTHEReigen
223*
224* =====================================================================
225 SUBROUTINE sstevx( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
226 $ M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
227*
228* -- LAPACK driver routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 CHARACTER JOBZ, RANGE
234 INTEGER IL, INFO, IU, LDZ, M, N
235 REAL ABSTOL, VL, VU
236* ..
237* .. Array Arguments ..
238 INTEGER IFAIL( * ), IWORK( * )
239 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 REAL ZERO, ONE
246 parameter( zero = 0.0e0, one = 1.0e0 )
247* ..
248* .. Local Scalars ..
249 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
250 CHARACTER ORDER
251 INTEGER I, IMAX, INDIBL, INDISP, INDIWO, INDWRK,
252 $ iscale, itmp1, j, jj, nsplit
253 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
254 $ tmp1, tnrm, vll, vuu
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 REAL SLAMCH, SLANST
259 EXTERNAL lsame, slamch, slanst
260* ..
261* .. External Subroutines ..
262 EXTERNAL scopy, sscal, sstebz, sstein, ssteqr, ssterf,
263 $ sswap, xerbla
264* ..
265* .. Intrinsic Functions ..
266 INTRINSIC max, min, sqrt
267* ..
268* .. Executable Statements ..
269*
270* Test the input parameters.
271*
272 wantz = lsame( jobz, 'V' )
273 alleig = lsame( range, 'A' )
274 valeig = lsame( range, 'V' )
275 indeig = lsame( range, 'I' )
276*
277 info = 0
278 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
279 info = -1
280 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
281 info = -2
282 ELSE IF( n.LT.0 ) THEN
283 info = -3
284 ELSE
285 IF( valeig ) THEN
286 IF( n.GT.0 .AND. vu.LE.vl )
287 $ info = -7
288 ELSE IF( indeig ) THEN
289 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
290 info = -8
291 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
292 info = -9
293 END IF
294 END IF
295 END IF
296 IF( info.EQ.0 ) THEN
297 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
298 $ info = -14
299 END IF
300*
301 IF( info.NE.0 ) THEN
302 CALL xerbla( 'SSTEVX', -info )
303 RETURN
304 END IF
305*
306* Quick return if possible
307*
308 m = 0
309 IF( n.EQ.0 )
310 $ RETURN
311*
312 IF( n.EQ.1 ) THEN
313 IF( alleig .OR. indeig ) THEN
314 m = 1
315 w( 1 ) = d( 1 )
316 ELSE
317 IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
318 m = 1
319 w( 1 ) = d( 1 )
320 END IF
321 END IF
322 IF( wantz )
323 $ z( 1, 1 ) = one
324 RETURN
325 END IF
326*
327* Get machine constants.
328*
329 safmin = slamch( 'Safe minimum' )
330 eps = slamch( 'Precision' )
331 smlnum = safmin / eps
332 bignum = one / smlnum
333 rmin = sqrt( smlnum )
334 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
335*
336* Scale matrix to allowable range, if necessary.
337*
338 iscale = 0
339 IF ( valeig ) THEN
340 vll = vl
341 vuu = vu
342 ELSE
343 vll = zero
344 vuu = zero
345 ENDIF
346 tnrm = slanst( 'M', n, d, e )
347 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
348 iscale = 1
349 sigma = rmin / tnrm
350 ELSE IF( tnrm.GT.rmax ) THEN
351 iscale = 1
352 sigma = rmax / tnrm
353 END IF
354 IF( iscale.EQ.1 ) THEN
355 CALL sscal( n, sigma, d, 1 )
356 CALL sscal( n-1, sigma, e( 1 ), 1 )
357 IF( valeig ) THEN
358 vll = vl*sigma
359 vuu = vu*sigma
360 END IF
361 END IF
362*
363* If all eigenvalues are desired and ABSTOL is less than zero, then
364* call SSTERF or SSTEQR. If this fails for some eigenvalue, then
365* try SSTEBZ.
366*
367 test = .false.
368 IF( indeig ) THEN
369 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
370 test = .true.
371 END IF
372 END IF
373 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
374 CALL scopy( n, d, 1, w, 1 )
375 CALL scopy( n-1, e( 1 ), 1, work( 1 ), 1 )
376 indwrk = n + 1
377 IF( .NOT.wantz ) THEN
378 CALL ssterf( n, w, work, info )
379 ELSE
380 CALL ssteqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
381 IF( info.EQ.0 ) THEN
382 DO 10 i = 1, n
383 ifail( i ) = 0
384 10 CONTINUE
385 END IF
386 END IF
387 IF( info.EQ.0 ) THEN
388 m = n
389 GO TO 20
390 END IF
391 info = 0
392 END IF
393*
394* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
395*
396 IF( wantz ) THEN
397 order = 'B'
398 ELSE
399 order = 'E'
400 END IF
401 indwrk = 1
402 indibl = 1
403 indisp = indibl + n
404 indiwo = indisp + n
405 CALL sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
406 $ nsplit, w, iwork( indibl ), iwork( indisp ),
407 $ work( indwrk ), iwork( indiwo ), info )
408*
409 IF( wantz ) THEN
410 CALL sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
411 $ z, ldz, work( indwrk ), iwork( indiwo ), ifail,
412 $ info )
413 END IF
414*
415* If matrix was scaled, then rescale eigenvalues appropriately.
416*
417 20 CONTINUE
418 IF( iscale.EQ.1 ) THEN
419 IF( info.EQ.0 ) THEN
420 imax = m
421 ELSE
422 imax = info - 1
423 END IF
424 CALL sscal( imax, one / sigma, w, 1 )
425 END IF
426*
427* If eigenvalues are not in order, then sort them, along with
428* eigenvectors.
429*
430 IF( wantz ) THEN
431 DO 40 j = 1, m - 1
432 i = 0
433 tmp1 = w( j )
434 DO 30 jj = j + 1, m
435 IF( w( jj ).LT.tmp1 ) THEN
436 i = jj
437 tmp1 = w( jj )
438 END IF
439 30 CONTINUE
440*
441 IF( i.NE.0 ) THEN
442 itmp1 = iwork( indibl+i-1 )
443 w( i ) = w( j )
444 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
445 w( j ) = tmp1
446 iwork( indibl+j-1 ) = itmp1
447 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
448 IF( info.NE.0 ) THEN
449 itmp1 = ifail( i )
450 ifail( i ) = ifail( j )
451 ifail( j ) = itmp1
452 END IF
453 END IF
454 40 CONTINUE
455 END IF
456*
457 RETURN
458*
459* End of SSTEVX
460*
461 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:273
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:174
subroutine sstevx(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: sstevx.f:227
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79