LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \endverbatim
93 *>
94 *> \param[in] VU
95 *> \verbatim
96 *> VU is REAL
97 *> If RANGE='V', the lower and upper bounds of the interval to
98 *> be searched for eigenvalues. VL < VU.
99 *> Not referenced if RANGE = 'A' or 'I'.
100 *> \endverbatim
101 *>
102 *> \param[in] IL
103 *> \verbatim
104 *> IL is INTEGER
105 *> \endverbatim
106 *>
107 *> \param[in] IU
108 *> \verbatim
109 *> IU is INTEGER
110 *> If RANGE='I', the indices (in ascending order) of the
111 *> smallest and largest eigenvalues to be returned.
112 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
113 *> Not referenced if RANGE = 'A' or 'V'.
114 *> \endverbatim
115 *>
116 *> \param[in] ABSTOL
117 *> \verbatim
118 *> ABSTOL is REAL
119 *> The absolute error tolerance for the eigenvalues.
120 *> An approximate eigenvalue is accepted as converged
121 *> when it is determined to lie in an interval [a,b]
122 *> of width less than or equal to
123 *>
124 *> ABSTOL + EPS * max( |a|,|b| ) ,
125 *>
126 *> where EPS is the machine precision. If ABSTOL is less
127 *> than or equal to zero, then EPS*|T| will be used in
128 *> its place, where |T| is the 1-norm of the tridiagonal
129 *> matrix.
130 *>
131 *> Eigenvalues will be computed most accurately when ABSTOL is
132 *> set to twice the underflow threshold 2*SLAMCH('S'), not zero.
133 *> If this routine returns with INFO>0, indicating that some
134 *> eigenvectors did not converge, try setting ABSTOL to
135 *> 2*SLAMCH('S').
136 *>
137 *> See "Computing Small Singular Values of Bidiagonal Matrices
138 *> with Guaranteed High Relative Accuracy," by Demmel and
139 *> Kahan, LAPACK Working Note #3.
140 *> \endverbatim
141 *>
142 *> \param[out] M
143 *> \verbatim
144 *> M is INTEGER
145 *> The total number of eigenvalues found. 0 <= M <= N.
146 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
147 *> \endverbatim
148 *>
149 *> \param[out] W
150 *> \verbatim
151 *> W is REAL array, dimension (N)
152 *> The first M elements contain the selected eigenvalues in
153 *> ascending order.
154 *> \endverbatim
155 *>
156 *> \param[out] Z
157 *> \verbatim
158 *> Z is REAL array, dimension (LDZ, max(1,M) )
159 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
160 *> contain the orthonormal eigenvectors of the matrix A
161 *> corresponding to the selected eigenvalues, with the i-th
162 *> column of Z holding the eigenvector associated with W(i).
163 *> If an eigenvector fails to converge (INFO > 0), then that
164 *> column of Z contains the latest approximation to the
165 *> eigenvector, and the index of the eigenvector is returned
166 *> in IFAIL. If JOBZ = 'N', then Z is not referenced.
167 *> Note: the user must ensure that at least max(1,M) columns are
168 *> supplied in the array Z; if RANGE = 'V', the exact value of M
169 *> is not known in advance and an upper bound must be used.
170 *> \endverbatim
171 *>
172 *> \param[in] LDZ
173 *> \verbatim
174 *> LDZ is INTEGER
175 *> The leading dimension of the array Z. LDZ >= 1, and if
176 *> JOBZ = 'V', LDZ >= max(1,N).
177 *> \endverbatim
178 *>
179 *> \param[out] WORK
180 *> \verbatim
181 *> WORK is REAL array, dimension (5*N)
182 *> \endverbatim
183 *>
184 *> \param[out] IWORK
185 *> \verbatim
186 *> IWORK is INTEGER array, dimension (5*N)
187 *> \endverbatim
188 *>
189 *> \param[out] IFAIL
190 *> \verbatim
191 *> IFAIL is INTEGER array, dimension (N)
192 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
193 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
194 *> indices of the eigenvectors that failed to converge.
195 *> If JOBZ = 'N', then IFAIL is not referenced.
196 *> \endverbatim
197 *>
198 *> \param[out] INFO
199 *> \verbatim
200 *> INFO is INTEGER
201 *> = 0: successful exit
202 *> < 0: if INFO = -i, the i-th argument had an illegal value
203 *> > 0: if INFO = i, then i eigenvectors failed to converge.
204 *> Their indices are stored in array IFAIL.
205 *> \endverbatim
206 *
207 * Authors:
208 * ========
209 *
210 *> \author Univ. of Tennessee
211 *> \author Univ. of California Berkeley
212 *> \author Univ. of Colorado Denver
213 *> \author NAG Ltd.
214 *
215 *> \date November 2011
216 *
217 *> \ingroup realOTHEReigen
218 *
219 * =====================================================================
220  SUBROUTINE sstevx( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
221  $ m, w, z, ldz, work, iwork, ifail, info )
222 *
223 * -- LAPACK driver routine (version 3.4.0) --
224 * -- LAPACK is a software package provided by Univ. of Tennessee, --
225 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226 * November 2011
227 *
228 * .. Scalar Arguments ..
229  CHARACTER jobz, range
230  INTEGER il, info, iu, ldz, m, n
231  REAL abstol, vl, vu
232 * ..
233 * .. Array Arguments ..
234  INTEGER ifail( * ), iwork( * )
235  REAL d( * ), e( * ), w( * ), work( * ), z( ldz, * )
236 * ..
237 *
238 * =====================================================================
239 *
240 * .. Parameters ..
241  REAL zero, one
242  parameter( zero = 0.0e0, one = 1.0e0 )
243 * ..
244 * .. Local Scalars ..
245  LOGICAL alleig, indeig, test, valeig, wantz
246  CHARACTER order
247  INTEGER i, imax, indibl, indisp, indiwo, indwrk,
248  $ iscale, itmp1, j, jj, nsplit
249  REAL bignum, eps, rmax, rmin, safmin, sigma, smlnum,
250  $ tmp1, tnrm, vll, vuu
251 * ..
252 * .. External Functions ..
253  LOGICAL lsame
254  REAL slamch, slanst
255  EXTERNAL lsame, slamch, slanst
256 * ..
257 * .. External Subroutines ..
258  EXTERNAL scopy, sscal, sstebz, sstein, ssteqr, ssterf,
259  $ sswap, xerbla
260 * ..
261 * .. Intrinsic Functions ..
262  INTRINSIC max, min, sqrt
263 * ..
264 * .. Executable Statements ..
265 *
266 * Test the input parameters.
267 *
268  wantz = lsame( jobz, 'V' )
269  alleig = lsame( range, 'A' )
270  valeig = lsame( range, 'V' )
271  indeig = lsame( range, 'I' )
272 *
273  info = 0
274  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
275  info = -1
276  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
277  info = -2
278  ELSE IF( n.LT.0 ) THEN
279  info = -3
280  ELSE
281  IF( valeig ) THEN
282  IF( n.GT.0 .AND. vu.LE.vl )
283  $ info = -7
284  ELSE IF( indeig ) THEN
285  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
286  info = -8
287  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
288  info = -9
289  END IF
290  END IF
291  END IF
292  IF( info.EQ.0 ) THEN
293  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
294  $ info = -14
295  END IF
296 *
297  IF( info.NE.0 ) THEN
298  CALL xerbla( 'SSTEVX', -info )
299  return
300  END IF
301 *
302 * Quick return if possible
303 *
304  m = 0
305  IF( n.EQ.0 )
306  $ return
307 *
308  IF( n.EQ.1 ) THEN
309  IF( alleig .OR. indeig ) THEN
310  m = 1
311  w( 1 ) = d( 1 )
312  ELSE
313  IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
314  m = 1
315  w( 1 ) = d( 1 )
316  END IF
317  END IF
318  IF( wantz )
319  $ z( 1, 1 ) = one
320  return
321  END IF
322 *
323 * Get machine constants.
324 *
325  safmin = slamch( 'Safe minimum' )
326  eps = slamch( 'Precision' )
327  smlnum = safmin / eps
328  bignum = one / smlnum
329  rmin = sqrt( smlnum )
330  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
331 *
332 * Scale matrix to allowable range, if necessary.
333 *
334  iscale = 0
335  IF ( valeig ) THEN
336  vll = vl
337  vuu = vu
338  ELSE
339  vll = zero
340  vuu = zero
341  ENDIF
342  tnrm = slanst( 'M', n, d, e )
343  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
344  iscale = 1
345  sigma = rmin / tnrm
346  ELSE IF( tnrm.GT.rmax ) THEN
347  iscale = 1
348  sigma = rmax / tnrm
349  END IF
350  IF( iscale.EQ.1 ) THEN
351  CALL sscal( n, sigma, d, 1 )
352  CALL sscal( n-1, sigma, e( 1 ), 1 )
353  IF( valeig ) THEN
354  vll = vl*sigma
355  vuu = vu*sigma
356  END IF
357  END IF
358 *
359 * If all eigenvalues are desired and ABSTOL is less than zero, then
360 * call SSTERF or SSTEQR. If this fails for some eigenvalue, then
361 * try SSTEBZ.
362 *
363  test = .false.
364  IF( indeig ) THEN
365  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
366  test = .true.
367  END IF
368  END IF
369  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
370  CALL scopy( n, d, 1, w, 1 )
371  CALL scopy( n-1, e( 1 ), 1, work( 1 ), 1 )
372  indwrk = n + 1
373  IF( .NOT.wantz ) THEN
374  CALL ssterf( n, w, work, info )
375  ELSE
376  CALL ssteqr( 'I', n, w, work, z, ldz, work( indwrk ), info )
377  IF( info.EQ.0 ) THEN
378  DO 10 i = 1, n
379  ifail( i ) = 0
380  10 continue
381  END IF
382  END IF
383  IF( info.EQ.0 ) THEN
384  m = n
385  go to 20
386  END IF
387  info = 0
388  END IF
389 *
390 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
391 *
392  IF( wantz ) THEN
393  order = 'B'
394  ELSE
395  order = 'E'
396  END IF
397  indwrk = 1
398  indibl = 1
399  indisp = indibl + n
400  indiwo = indisp + n
401  CALL sstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
402  $ nsplit, w, iwork( indibl ), iwork( indisp ),
403  $ work( indwrk ), iwork( indiwo ), info )
404 *
405  IF( wantz ) THEN
406  CALL sstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
407  $ z, ldz, work( indwrk ), iwork( indiwo ), ifail,
408  $ info )
409  END IF
410 *
411 * If matrix was scaled, then rescale eigenvalues appropriately.
412 *
413  20 continue
414  IF( iscale.EQ.1 ) THEN
415  IF( info.EQ.0 ) THEN
416  imax = m
417  ELSE
418  imax = info - 1
419  END IF
420  CALL sscal( imax, one / sigma, w, 1 )
421  END IF
422 *
423 * If eigenvalues are not in order, then sort them, along with
424 * eigenvectors.
425 *
426  IF( wantz ) THEN
427  DO 40 j = 1, m - 1
428  i = 0
429  tmp1 = w( j )
430  DO 30 jj = j + 1, m
431  IF( w( jj ).LT.tmp1 ) THEN
432  i = jj
433  tmp1 = w( jj )
434  END IF
435  30 continue
436 *
437  IF( i.NE.0 ) THEN
438  itmp1 = iwork( indibl+i-1 )
439  w( i ) = w( j )
440  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
441  w( j ) = tmp1
442  iwork( indibl+j-1 ) = itmp1
443  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
444  IF( info.NE.0 ) THEN
445  itmp1 = ifail( i )
446  ifail( i ) = ifail( j )
447  ifail( j ) = itmp1
448  END IF
449  END IF
450  40 continue
451  END IF
452 *
453  return
454 *
455 * End of SSTEVX
456 *
457  END