LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
sstemr.f
Go to the documentation of this file.
1 *> \brief \b SSTEMR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSTEMR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstemr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstemr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstemr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
22 * M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
23 * IWORK, LIWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE
27 * LOGICAL TRYRAC
28 * INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
29 * REAL VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * ), IWORK( * )
33 * REAL D( * ), E( * ), W( * ), WORK( * )
34 * REAL Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> SSTEMR computes selected eigenvalues and, optionally, eigenvectors
44 *> of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
45 *> a well defined set of pairwise different real eigenvalues, the corresponding
46 *> real eigenvectors are pairwise orthogonal.
47 *>
48 *> The spectrum may be computed either completely or partially by specifying
49 *> either an interval (VL,VU] or a range of indices IL:IU for the desired
50 *> eigenvalues.
51 *>
52 *> Depending on the number of desired eigenvalues, these are computed either
53 *> by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
54 *> computed by the use of various suitable L D L^T factorizations near clusters
55 *> of close eigenvalues (referred to as RRRs, Relatively Robust
56 *> Representations). An informal sketch of the algorithm follows.
57 *>
58 *> For each unreduced block (submatrix) of T,
59 *> (a) Compute T - sigma I = L D L^T, so that L and D
60 *> define all the wanted eigenvalues to high relative accuracy.
61 *> This means that small relative changes in the entries of D and L
62 *> cause only small relative changes in the eigenvalues and
63 *> eigenvectors. The standard (unfactored) representation of the
64 *> tridiagonal matrix T does not have this property in general.
65 *> (b) Compute the eigenvalues to suitable accuracy.
66 *> If the eigenvectors are desired, the algorithm attains full
67 *> accuracy of the computed eigenvalues only right before
68 *> the corresponding vectors have to be computed, see steps c) and d).
69 *> (c) For each cluster of close eigenvalues, select a new
70 *> shift close to the cluster, find a new factorization, and refine
71 *> the shifted eigenvalues to suitable accuracy.
72 *> (d) For each eigenvalue with a large enough relative separation compute
73 *> the corresponding eigenvector by forming a rank revealing twisted
74 *> factorization. Go back to (c) for any clusters that remain.
75 *>
76 *> For more details, see:
77 *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
78 *> to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
79 *> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
80 *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
81 *> Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
82 *> 2004. Also LAPACK Working Note 154.
83 *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
84 *> tridiagonal eigenvalue/eigenvector problem",
85 *> Computer Science Division Technical Report No. UCB/CSD-97-971,
86 *> UC Berkeley, May 1997.
87 *>
88 *> Further Details
89 *> 1.SSTEMR works only on machines which follow IEEE-754
90 *> floating-point standard in their handling of infinities and NaNs.
91 *> This permits the use of efficient inner loops avoiding a check for
92 *> zero divisors.
93 *> \endverbatim
94 *
95 * Arguments:
96 * ==========
97 *
98 *> \param[in] JOBZ
99 *> \verbatim
100 *> JOBZ is CHARACTER*1
101 *> = 'N': Compute eigenvalues only;
102 *> = 'V': Compute eigenvalues and eigenvectors.
103 *> \endverbatim
104 *>
105 *> \param[in] RANGE
106 *> \verbatim
107 *> RANGE is CHARACTER*1
108 *> = 'A': all eigenvalues will be found.
109 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
110 *> will be found.
111 *> = 'I': the IL-th through IU-th eigenvalues will be found.
112 *> \endverbatim
113 *>
114 *> \param[in] N
115 *> \verbatim
116 *> N is INTEGER
117 *> The order of the matrix. N >= 0.
118 *> \endverbatim
119 *>
120 *> \param[in,out] D
121 *> \verbatim
122 *> D is REAL array, dimension (N)
123 *> On entry, the N diagonal elements of the tridiagonal matrix
124 *> T. On exit, D is overwritten.
125 *> \endverbatim
126 *>
127 *> \param[in,out] E
128 *> \verbatim
129 *> E is REAL array, dimension (N)
130 *> On entry, the (N-1) subdiagonal elements of the tridiagonal
131 *> matrix T in elements 1 to N-1 of E. E(N) need not be set on
132 *> input, but is used internally as workspace.
133 *> On exit, E is overwritten.
134 *> \endverbatim
135 *>
136 *> \param[in] VL
137 *> \verbatim
138 *> VL is REAL
139 *> \endverbatim
140 *>
141 *> \param[in] VU
142 *> \verbatim
143 *> VU is REAL
144 *>
145 *> If RANGE='V', the lower and upper bounds of the interval to
146 *> be searched for eigenvalues. VL < VU.
147 *> Not referenced if RANGE = 'A' or 'I'.
148 *> \endverbatim
149 *>
150 *> \param[in] IL
151 *> \verbatim
152 *> IL is INTEGER
153 *> \endverbatim
154 *>
155 *> \param[in] IU
156 *> \verbatim
157 *> IU is INTEGER
158 *>
159 *> If RANGE='I', the indices (in ascending order) of the
160 *> smallest and largest eigenvalues to be returned.
161 *> 1 <= IL <= IU <= N, if N > 0.
162 *> Not referenced if RANGE = 'A' or 'V'.
163 *> \endverbatim
164 *>
165 *> \param[out] M
166 *> \verbatim
167 *> M is INTEGER
168 *> The total number of eigenvalues found. 0 <= M <= N.
169 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
170 *> \endverbatim
171 *>
172 *> \param[out] W
173 *> \verbatim
174 *> W is REAL array, dimension (N)
175 *> The first M elements contain the selected eigenvalues in
176 *> ascending order.
177 *> \endverbatim
178 *>
179 *> \param[out] Z
180 *> \verbatim
181 *> Z is REAL array, dimension (LDZ, max(1,M) )
182 *> If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
183 *> contain the orthonormal eigenvectors of the matrix T
184 *> corresponding to the selected eigenvalues, with the i-th
185 *> column of Z holding the eigenvector associated with W(i).
186 *> If JOBZ = 'N', then Z is not referenced.
187 *> Note: the user must ensure that at least max(1,M) columns are
188 *> supplied in the array Z; if RANGE = 'V', the exact value of M
189 *> is not known in advance and can be computed with a workspace
190 *> query by setting NZC = -1, see below.
191 *> \endverbatim
192 *>
193 *> \param[in] LDZ
194 *> \verbatim
195 *> LDZ is INTEGER
196 *> The leading dimension of the array Z. LDZ >= 1, and if
197 *> JOBZ = 'V', then LDZ >= max(1,N).
198 *> \endverbatim
199 *>
200 *> \param[in] NZC
201 *> \verbatim
202 *> NZC is INTEGER
203 *> The number of eigenvectors to be held in the array Z.
204 *> If RANGE = 'A', then NZC >= max(1,N).
205 *> If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
206 *> If RANGE = 'I', then NZC >= IU-IL+1.
207 *> If NZC = -1, then a workspace query is assumed; the
208 *> routine calculates the number of columns of the array Z that
209 *> are needed to hold the eigenvectors.
210 *> This value is returned as the first entry of the Z array, and
211 *> no error message related to NZC is issued by XERBLA.
212 *> \endverbatim
213 *>
214 *> \param[out] ISUPPZ
215 *> \verbatim
216 *> ISUPPZ is INTEGER ARRAY, dimension ( 2*max(1,M) )
217 *> The support of the eigenvectors in Z, i.e., the indices
218 *> indicating the nonzero elements in Z. The i-th computed eigenvector
219 *> is nonzero only in elements ISUPPZ( 2*i-1 ) through
220 *> ISUPPZ( 2*i ). This is relevant in the case when the matrix
221 *> is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
222 *> \endverbatim
223 *>
224 *> \param[in,out] TRYRAC
225 *> \verbatim
226 *> TRYRAC is LOGICAL
227 *> If TRYRAC.EQ..TRUE., indicates that the code should check whether
228 *> the tridiagonal matrix defines its eigenvalues to high relative
229 *> accuracy. If so, the code uses relative-accuracy preserving
230 *> algorithms that might be (a bit) slower depending on the matrix.
231 *> If the matrix does not define its eigenvalues to high relative
232 *> accuracy, the code can uses possibly faster algorithms.
233 *> If TRYRAC.EQ..FALSE., the code is not required to guarantee
234 *> relatively accurate eigenvalues and can use the fastest possible
235 *> techniques.
236 *> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
237 *> does not define its eigenvalues to high relative accuracy.
238 *> \endverbatim
239 *>
240 *> \param[out] WORK
241 *> \verbatim
242 *> WORK is REAL array, dimension (LWORK)
243 *> On exit, if INFO = 0, WORK(1) returns the optimal
244 *> (and minimal) LWORK.
245 *> \endverbatim
246 *>
247 *> \param[in] LWORK
248 *> \verbatim
249 *> LWORK is INTEGER
250 *> The dimension of the array WORK. LWORK >= max(1,18*N)
251 *> if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
252 *> If LWORK = -1, then a workspace query is assumed; the routine
253 *> only calculates the optimal size of the WORK array, returns
254 *> this value as the first entry of the WORK array, and no error
255 *> message related to LWORK is issued by XERBLA.
256 *> \endverbatim
257 *>
258 *> \param[out] IWORK
259 *> \verbatim
260 *> IWORK is INTEGER array, dimension (LIWORK)
261 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
262 *> \endverbatim
263 *>
264 *> \param[in] LIWORK
265 *> \verbatim
266 *> LIWORK is INTEGER
267 *> The dimension of the array IWORK. LIWORK >= max(1,10*N)
268 *> if the eigenvectors are desired, and LIWORK >= max(1,8*N)
269 *> if only the eigenvalues are to be computed.
270 *> If LIWORK = -1, then a workspace query is assumed; the
271 *> routine only calculates the optimal size of the IWORK array,
272 *> returns this value as the first entry of the IWORK array, and
273 *> no error message related to LIWORK is issued by XERBLA.
274 *> \endverbatim
275 *>
276 *> \param[out] INFO
277 *> \verbatim
278 *> INFO is INTEGER
279 *> On exit, INFO
280 *> = 0: successful exit
281 *> < 0: if INFO = -i, the i-th argument had an illegal value
282 *> > 0: if INFO = 1X, internal error in SLARRE,
283 *> if INFO = 2X, internal error in SLARRV.
284 *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
285 *> the nonzero error code returned by SLARRE or
286 *> SLARRV, respectively.
287 *> \endverbatim
288 *
289 * Authors:
290 * ========
291 *
292 *> \author Univ. of Tennessee
293 *> \author Univ. of California Berkeley
294 *> \author Univ. of Colorado Denver
295 *> \author NAG Ltd.
296 *
297 *> \date November 2013
298 *
299 *> \ingroup realOTHERcomputational
300 *
301 *> \par Contributors:
302 * ==================
303 *>
304 *> Beresford Parlett, University of California, Berkeley, USA \n
305 *> Jim Demmel, University of California, Berkeley, USA \n
306 *> Inderjit Dhillon, University of Texas, Austin, USA \n
307 *> Osni Marques, LBNL/NERSC, USA \n
308 *> Christof Voemel, University of California, Berkeley, USA
309 *
310 * =====================================================================
311  SUBROUTINE sstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
312  $ m, w, z, ldz, nzc, isuppz, tryrac, work, lwork,
313  $ iwork, liwork, info )
314 *
315 * -- LAPACK computational routine (version 3.5.0) --
316 * -- LAPACK is a software package provided by Univ. of Tennessee, --
317 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
318 * November 2013
319 *
320 * .. Scalar Arguments ..
321  CHARACTER jobz, range
322  LOGICAL tryrac
323  INTEGER il, info, iu, ldz, nzc, liwork, lwork, m, n
324  REAL vl, vu
325 * ..
326 * .. Array Arguments ..
327  INTEGER isuppz( * ), iwork( * )
328  REAL d( * ), e( * ), w( * ), work( * )
329  REAL z( ldz, * )
330 * ..
331 *
332 * =====================================================================
333 *
334 * .. Parameters ..
335  REAL zero, one, four, minrgp
336  parameter( zero = 0.0e0, one = 1.0e0,
337  $ four = 4.0e0,
338  $ minrgp = 3.0e-3 )
339 * ..
340 * .. Local Scalars ..
341  LOGICAL alleig, indeig, lquery, valeig, wantz, zquery
342  INTEGER i, ibegin, iend, ifirst, iil, iindbl, iindw,
343  $ iindwk, iinfo, iinspl, iiu, ilast, in, indd,
344  $ inde2, inderr, indgp, indgrs, indwrk, itmp,
345  $ itmp2, j, jblk, jj, liwmin, lwmin, nsplit,
346  $ nzcmin, offset, wbegin, wend
347  REAL bignum, cs, eps, pivmin, r1, r2, rmax, rmin,
348  $ rtol1, rtol2, safmin, scale, smlnum, sn,
349  $ thresh, tmp, tnrm, wl, wu
350 * ..
351 * ..
352 * .. External Functions ..
353  LOGICAL lsame
354  REAL slamch, slanst
355  EXTERNAL lsame, slamch, slanst
356 * ..
357 * .. External Subroutines ..
358  EXTERNAL scopy, slae2, slaev2, slarrc, slarre, slarrj,
360 * ..
361 * .. Intrinsic Functions ..
362  INTRINSIC max, min, sqrt
363 * ..
364 * .. Executable Statements ..
365 *
366 * Test the input parameters.
367 *
368  wantz = lsame( jobz, 'V' )
369  alleig = lsame( range, 'A' )
370  valeig = lsame( range, 'V' )
371  indeig = lsame( range, 'I' )
372 *
373  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
374  zquery = ( nzc.EQ.-1 )
375 
376 * SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
377 * In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
378 * Furthermore, SLARRV needs WORK of size 12*N, IWORK of size 7*N.
379  IF( wantz ) THEN
380  lwmin = 18*n
381  liwmin = 10*n
382  ELSE
383 * need less workspace if only the eigenvalues are wanted
384  lwmin = 12*n
385  liwmin = 8*n
386  ENDIF
387 
388  wl = zero
389  wu = zero
390  iil = 0
391  iiu = 0
392  nsplit = 0
393 
394  IF( valeig ) THEN
395 * We do not reference VL, VU in the cases RANGE = 'I','A'
396 * The interval (WL, WU] contains all the wanted eigenvalues.
397 * It is either given by the user or computed in SLARRE.
398  wl = vl
399  wu = vu
400  ELSEIF( indeig ) THEN
401 * We do not reference IL, IU in the cases RANGE = 'V','A'
402  iil = il
403  iiu = iu
404  ENDIF
405 *
406  info = 0
407  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
408  info = -1
409  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
410  info = -2
411  ELSE IF( n.LT.0 ) THEN
412  info = -3
413  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
414  info = -7
415  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
416  info = -8
417  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
418  info = -9
419  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
420  info = -13
421  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
422  info = -17
423  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
424  info = -19
425  END IF
426 *
427 * Get machine constants.
428 *
429  safmin = slamch( 'Safe minimum' )
430  eps = slamch( 'Precision' )
431  smlnum = safmin / eps
432  bignum = one / smlnum
433  rmin = sqrt( smlnum )
434  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
435 *
436  IF( info.EQ.0 ) THEN
437  work( 1 ) = lwmin
438  iwork( 1 ) = liwmin
439 *
440  IF( wantz .AND. alleig ) THEN
441  nzcmin = n
442  ELSE IF( wantz .AND. valeig ) THEN
443  CALL slarrc( 'T', n, vl, vu, d, e, safmin,
444  $ nzcmin, itmp, itmp2, info )
445  ELSE IF( wantz .AND. indeig ) THEN
446  nzcmin = iiu-iil+1
447  ELSE
448 * WANTZ .EQ. FALSE.
449  nzcmin = 0
450  ENDIF
451  IF( zquery .AND. info.EQ.0 ) THEN
452  z( 1,1 ) = nzcmin
453  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
454  info = -14
455  END IF
456  END IF
457 
458  IF( info.NE.0 ) THEN
459 *
460  CALL xerbla( 'SSTEMR', -info )
461 *
462  RETURN
463  ELSE IF( lquery .OR. zquery ) THEN
464  RETURN
465  END IF
466 *
467 * Handle N = 0, 1, and 2 cases immediately
468 *
469  m = 0
470  IF( n.EQ.0 )
471  $ RETURN
472 *
473  IF( n.EQ.1 ) THEN
474  IF( alleig .OR. indeig ) THEN
475  m = 1
476  w( 1 ) = d( 1 )
477  ELSE
478  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
479  m = 1
480  w( 1 ) = d( 1 )
481  END IF
482  END IF
483  IF( wantz.AND.(.NOT.zquery) ) THEN
484  z( 1, 1 ) = one
485  isuppz(1) = 1
486  isuppz(2) = 1
487  END IF
488  RETURN
489  END IF
490 *
491  IF( n.EQ.2 ) THEN
492  IF( .NOT.wantz ) THEN
493  CALL slae2( d(1), e(1), d(2), r1, r2 )
494  ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
495  CALL slaev2( d(1), e(1), d(2), r1, r2, cs, sn )
496  END IF
497  IF( alleig.OR.
498  $ (valeig.AND.(r2.GT.wl).AND.
499  $ (r2.LE.wu)).OR.
500  $ (indeig.AND.(iil.EQ.1)) ) THEN
501  m = m+1
502  w( m ) = r2
503  IF( wantz.AND.(.NOT.zquery) ) THEN
504  z( 1, m ) = -sn
505  z( 2, m ) = cs
506 * Note: At most one of SN and CS can be zero.
507  IF (sn.NE.zero) THEN
508  IF (cs.NE.zero) THEN
509  isuppz(2*m-1) = 1
510  isuppz(2*m) = 2
511  ELSE
512  isuppz(2*m-1) = 1
513  isuppz(2*m) = 1
514  END IF
515  ELSE
516  isuppz(2*m-1) = 2
517  isuppz(2*m) = 2
518  END IF
519  ENDIF
520  ENDIF
521  IF( alleig.OR.
522  $ (valeig.AND.(r1.GT.wl).AND.
523  $ (r1.LE.wu)).OR.
524  $ (indeig.AND.(iiu.EQ.2)) ) THEN
525  m = m+1
526  w( m ) = r1
527  IF( wantz.AND.(.NOT.zquery) ) THEN
528  z( 1, m ) = cs
529  z( 2, m ) = sn
530 * Note: At most one of SN and CS can be zero.
531  IF (sn.NE.zero) THEN
532  IF (cs.NE.zero) THEN
533  isuppz(2*m-1) = 1
534  isuppz(2*m) = 2
535  ELSE
536  isuppz(2*m-1) = 1
537  isuppz(2*m) = 1
538  END IF
539  ELSE
540  isuppz(2*m-1) = 2
541  isuppz(2*m) = 2
542  END IF
543  ENDIF
544  ENDIF
545  ELSE
546 
547 * Continue with general N
548 
549  indgrs = 1
550  inderr = 2*n + 1
551  indgp = 3*n + 1
552  indd = 4*n + 1
553  inde2 = 5*n + 1
554  indwrk = 6*n + 1
555 *
556  iinspl = 1
557  iindbl = n + 1
558  iindw = 2*n + 1
559  iindwk = 3*n + 1
560 *
561 * Scale matrix to allowable range, if necessary.
562 * The allowable range is related to the PIVMIN parameter; see the
563 * comments in SLARRD. The preference for scaling small values
564 * up is heuristic; we expect users' matrices not to be close to the
565 * RMAX threshold.
566 *
567  scale = one
568  tnrm = slanst( 'M', n, d, e )
569  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
570  scale = rmin / tnrm
571  ELSE IF( tnrm.GT.rmax ) THEN
572  scale = rmax / tnrm
573  END IF
574  IF( scale.NE.one ) THEN
575  CALL sscal( n, scale, d, 1 )
576  CALL sscal( n-1, scale, e, 1 )
577  tnrm = tnrm*scale
578  IF( valeig ) THEN
579 * If eigenvalues in interval have to be found,
580 * scale (WL, WU] accordingly
581  wl = wl*scale
582  wu = wu*scale
583  ENDIF
584  END IF
585 *
586 * Compute the desired eigenvalues of the tridiagonal after splitting
587 * into smaller subblocks if the corresponding off-diagonal elements
588 * are small
589 * THRESH is the splitting parameter for SLARRE
590 * A negative THRESH forces the old splitting criterion based on the
591 * size of the off-diagonal. A positive THRESH switches to splitting
592 * which preserves relative accuracy.
593 *
594  IF( tryrac ) THEN
595 * Test whether the matrix warrants the more expensive relative approach.
596  CALL slarrr( n, d, e, iinfo )
597  ELSE
598 * The user does not care about relative accurately eigenvalues
599  iinfo = -1
600  ENDIF
601 * Set the splitting criterion
602  IF (iinfo.EQ.0) THEN
603  thresh = eps
604  ELSE
605  thresh = -eps
606 * relative accuracy is desired but T does not guarantee it
607  tryrac = .false.
608  ENDIF
609 *
610  IF( tryrac ) THEN
611 * Copy original diagonal, needed to guarantee relative accuracy
612  CALL scopy(n,d,1,work(indd),1)
613  ENDIF
614 * Store the squares of the offdiagonal values of T
615  DO 5 j = 1, n-1
616  work( inde2+j-1 ) = e(j)**2
617  5 CONTINUE
618 
619 * Set the tolerance parameters for bisection
620  IF( .NOT.wantz ) THEN
621 * SLARRE computes the eigenvalues to full precision.
622  rtol1 = four * eps
623  rtol2 = four * eps
624  ELSE
625 * SLARRE computes the eigenvalues to less than full precision.
626 * SLARRV will refine the eigenvalue approximations, and we can
627 * need less accurate initial bisection in SLARRE.
628 * Note: these settings do only affect the subset case and SLARRE
629  rtol1 = max( sqrt(eps)*5.0e-2, four * eps )
630  rtol2 = max( sqrt(eps)*5.0e-3, four * eps )
631  ENDIF
632  CALL slarre( range, n, wl, wu, iil, iiu, d, e,
633  $ work(inde2), rtol1, rtol2, thresh, nsplit,
634  $ iwork( iinspl ), m, w, work( inderr ),
635  $ work( indgp ), iwork( iindbl ),
636  $ iwork( iindw ), work( indgrs ), pivmin,
637  $ work( indwrk ), iwork( iindwk ), iinfo )
638  IF( iinfo.NE.0 ) THEN
639  info = 10 + abs( iinfo )
640  RETURN
641  END IF
642 * Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
643 * part of the spectrum. All desired eigenvalues are contained in
644 * (WL,WU]
645 
646 
647  IF( wantz ) THEN
648 *
649 * Compute the desired eigenvectors corresponding to the computed
650 * eigenvalues
651 *
652  CALL slarrv( n, wl, wu, d, e,
653  $ pivmin, iwork( iinspl ), m,
654  $ 1, m, minrgp, rtol1, rtol2,
655  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
656  $ iwork( iindw ), work( indgrs ), z, ldz,
657  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
658  IF( iinfo.NE.0 ) THEN
659  info = 20 + abs( iinfo )
660  RETURN
661  END IF
662  ELSE
663 * SLARRE computes eigenvalues of the (shifted) root representation
664 * SLARRV returns the eigenvalues of the unshifted matrix.
665 * However, if the eigenvectors are not desired by the user, we need
666 * to apply the corresponding shifts from SLARRE to obtain the
667 * eigenvalues of the original matrix.
668  DO 20 j = 1, m
669  itmp = iwork( iindbl+j-1 )
670  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
671  20 CONTINUE
672  END IF
673 *
674 
675  IF ( tryrac ) THEN
676 * Refine computed eigenvalues so that they are relatively accurate
677 * with respect to the original matrix T.
678  ibegin = 1
679  wbegin = 1
680  DO 39 jblk = 1, iwork( iindbl+m-1 )
681  iend = iwork( iinspl+jblk-1 )
682  in = iend - ibegin + 1
683  wend = wbegin - 1
684 * check if any eigenvalues have to be refined in this block
685  36 CONTINUE
686  IF( wend.LT.m ) THEN
687  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
688  wend = wend + 1
689  go to 36
690  END IF
691  END IF
692  IF( wend.LT.wbegin ) THEN
693  ibegin = iend + 1
694  go to 39
695  END IF
696 
697  offset = iwork(iindw+wbegin-1)-1
698  ifirst = iwork(iindw+wbegin-1)
699  ilast = iwork(iindw+wend-1)
700  rtol2 = four * eps
701  CALL slarrj( in,
702  $ work(indd+ibegin-1), work(inde2+ibegin-1),
703  $ ifirst, ilast, rtol2, offset, w(wbegin),
704  $ work( inderr+wbegin-1 ),
705  $ work( indwrk ), iwork( iindwk ), pivmin,
706  $ tnrm, iinfo )
707  ibegin = iend + 1
708  wbegin = wend + 1
709  39 CONTINUE
710  ENDIF
711 *
712 * If matrix was scaled, then rescale eigenvalues appropriately.
713 *
714  IF( scale.NE.one ) THEN
715  CALL sscal( m, one / scale, w, 1 )
716  END IF
717  END IF
718 *
719 * If eigenvalues are not in increasing order, then sort them,
720 * possibly along with eigenvectors.
721 *
722  IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
723  IF( .NOT. wantz ) THEN
724  CALL slasrt( 'I', m, w, iinfo )
725  IF( iinfo.NE.0 ) THEN
726  info = 3
727  RETURN
728  END IF
729  ELSE
730  DO 60 j = 1, m - 1
731  i = 0
732  tmp = w( j )
733  DO 50 jj = j + 1, m
734  IF( w( jj ).LT.tmp ) THEN
735  i = jj
736  tmp = w( jj )
737  END IF
738  50 CONTINUE
739  IF( i.NE.0 ) THEN
740  w( i ) = w( j )
741  w( j ) = tmp
742  IF( wantz ) THEN
743  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
744  itmp = isuppz( 2*i-1 )
745  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
746  isuppz( 2*j-1 ) = itmp
747  itmp = isuppz( 2*i )
748  isuppz( 2*i ) = isuppz( 2*j )
749  isuppz( 2*j ) = itmp
750  END IF
751  END IF
752  60 CONTINUE
753  END IF
754  ENDIF
755 *
756 *
757  work( 1 ) = lwmin
758  iwork( 1 ) = liwmin
759  RETURN
760 *
761 * End of SSTEMR
762 *
763  END
subroutine slarrr(N, D, E, INFO)
SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computa...
Definition: slarrr.f:95
LOGICAL function lsame(CA, CB)
LSAME
Definition: lsame.f:54
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:52
REAL function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
REAL function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: slanst.f:101
subroutine slae2(A, B, C, RT1, RT2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: slae2.f:103
subroutine slaev2(A, B, C, RT1, RT2, CS1, SN1)
SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition: slaev2.f:121
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:52
subroutine slarrv(N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, WORK, IWORK, INFO)
SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition: slarrv.f:280
subroutine slarre(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
Definition: slarre.f:295
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine slarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
Definition: slarrc.f:136
subroutine sstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEMR
Definition: sstemr.f:311
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:89
subroutine slarrj(N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO)
SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T...
Definition: slarrj.f:167
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:54