LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dstevx.f
Go to the documentation of this file.
1*> \brief <b> DSTEVX 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 DSTEVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstevx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstevx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstevx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSTEVX( 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* DOUBLE PRECISION ABSTOL, VL, VU
28* ..
29* .. Array Arguments ..
30* INTEGER IFAIL( * ), IWORK( * )
31* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> DSTEVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION
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*DLAMCH('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*DLAMCH('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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 stevx
223*
224* =====================================================================
225 SUBROUTINE dstevx( 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 DOUBLE PRECISION ABSTOL, VL, VU
236* ..
237* .. Array Arguments ..
238 INTEGER IFAIL( * ), IWORK( * )
239 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
240* ..
241*
242* =====================================================================
243*
244* .. Parameters ..
245 DOUBLE PRECISION ZERO, ONE
246 parameter( zero = 0.0d0, one = 1.0d0 )
247* ..
248* .. Local Scalars ..
249 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
250 CHARACTER ORDER
251 INTEGER I, IMAX, INDISP, INDIWO, INDWRK,
252 $ iscale, itmp1, j, jj, nsplit
253 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
254 $ tmp1, tnrm, vll, vuu
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 DOUBLE PRECISION DLAMCH, DLANST
259 EXTERNAL lsame, dlamch, dlanst
260* ..
261* .. External Subroutines ..
262 EXTERNAL dcopy, dscal, dstebz, dstein, dsteqr, dsterf,
263 $ dswap, 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( 'DSTEVX', -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 = dlamch( 'Safe minimum' )
330 eps = dlamch( '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 END IF
346 tnrm = dlanst( '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 dscal( n, sigma, d, 1 )
356 CALL dscal( 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 DSTERF or SSTEQR. If this fails for some eigenvalue, then
365* try DSTEBZ.
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 dcopy( n, d, 1, w, 1 )
375 CALL dcopy( n-1, e( 1 ), 1, work( 1 ), 1 )
376 indwrk = n + 1
377 IF( .NOT.wantz ) THEN
378 CALL dsterf( n, w, work, info )
379 ELSE
380 CALL dsteqr( '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 DSTEBZ 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 indisp = 1 + n
403 indiwo = indisp + n
404 CALL dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
405 $ nsplit, w, iwork( 1 ), iwork( indisp ),
406 $ work( indwrk ), iwork( indiwo ), info )
407*
408 IF( wantz ) THEN
409 CALL dstein( n, d, e, m, w, iwork( 1 ), iwork( indisp ),
410 $ z, ldz, work( indwrk ), iwork( indiwo ), ifail,
411 $ info )
412 END IF
413*
414* If matrix was scaled, then rescale eigenvalues appropriately.
415*
416 20 CONTINUE
417 IF( iscale.EQ.1 ) THEN
418 IF( info.EQ.0 ) THEN
419 imax = m
420 ELSE
421 imax = info - 1
422 END IF
423 CALL dscal( imax, one / sigma, w, 1 )
424 END IF
425*
426* If eigenvalues are not in order, then sort them, along with
427* eigenvectors.
428*
429 IF( wantz ) THEN
430 DO 40 j = 1, m - 1
431 i = 0
432 tmp1 = w( j )
433 DO 30 jj = j + 1, m
434 IF( w( jj ).LT.tmp1 ) THEN
435 i = jj
436 tmp1 = w( jj )
437 END IF
438 30 CONTINUE
439*
440 IF( i.NE.0 ) THEN
441 itmp1 = iwork( 1 + i-1 )
442 w( i ) = w( j )
443 iwork( 1 + i-1 ) = iwork( 1 + j-1 )
444 w( j ) = tmp1
445 iwork( 1 + j-1 ) = itmp1
446 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
447 IF( info.NE.0 ) THEN
448 itmp1 = ifail( i )
449 ifail( i ) = ifail( j )
450 ifail( j ) = itmp1
451 END IF
452 END IF
453 40 CONTINUE
454 END IF
455*
456 RETURN
457*
458* End of DSTEVX
459*
460 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:273
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:174
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine dstevx(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
DSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition dstevx.f:227
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82