LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dstevx.f
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, INDIBL, 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  indibl = 1
403  indisp = indibl + n
404  indiwo = indisp + n
405  CALL dstebz( 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 dstein( 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 dscal( 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 dswap( 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 DSTEVX
460 *
461  END
dstebz
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
dsterf
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
dstein
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:174
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
dstevx
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 matric...
Definition: dstevx.f:227
dsteqr
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131