LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*>
140*> If RANGE='V', the lower bound of the interval to
141*> be searched for eigenvalues. VL < VU.
142*> Not referenced if RANGE = 'A' or 'I'.
143*> \endverbatim
144*>
145*> \param[in] VU
146*> \verbatim
147*> VU is REAL
148*>
149*> If RANGE='V', the upper bound of the interval to
150*> be searched for eigenvalues. VL < VU.
151*> Not referenced if RANGE = 'A' or 'I'.
152*> \endverbatim
153*>
154*> \param[in] IL
155*> \verbatim
156*> IL is INTEGER
157*>
158*> If RANGE='I', the index of the
159*> smallest eigenvalue to be returned.
160*> 1 <= IL <= IU <= N, if N > 0.
161*> Not referenced if RANGE = 'A' or 'V'.
162*> \endverbatim
163*>
164*> \param[in] IU
165*> \verbatim
166*> IU is INTEGER
167*>
168*> If RANGE='I', the index of the
169*> largest eigenvalue to be returned.
170*> 1 <= IL <= IU <= N, if N > 0.
171*> Not referenced if RANGE = 'A' or 'V'.
172*> \endverbatim
173*>
174*> \param[out] M
175*> \verbatim
176*> M is INTEGER
177*> The total number of eigenvalues found. 0 <= M <= N.
178*> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
179*> \endverbatim
180*>
181*> \param[out] W
182*> \verbatim
183*> W is REAL array, dimension (N)
184*> The first M elements contain the selected eigenvalues in
185*> ascending order.
186*> \endverbatim
187*>
188*> \param[out] Z
189*> \verbatim
190*> Z is REAL array, dimension (LDZ, max(1,M) )
191*> If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
192*> contain the orthonormal eigenvectors of the matrix T
193*> corresponding to the selected eigenvalues, with the i-th
194*> column of Z holding the eigenvector associated with W(i).
195*> If JOBZ = 'N', then Z is not referenced.
196*> Note: the user must ensure that at least max(1,M) columns are
197*> supplied in the array Z; if RANGE = 'V', the exact value of M
198*> is not known in advance and can be computed with a workspace
199*> query by setting NZC = -1, see below.
200*> \endverbatim
201*>
202*> \param[in] LDZ
203*> \verbatim
204*> LDZ is INTEGER
205*> The leading dimension of the array Z. LDZ >= 1, and if
206*> JOBZ = 'V', then LDZ >= max(1,N).
207*> \endverbatim
208*>
209*> \param[in] NZC
210*> \verbatim
211*> NZC is INTEGER
212*> The number of eigenvectors to be held in the array Z.
213*> If RANGE = 'A', then NZC >= max(1,N).
214*> If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
215*> If RANGE = 'I', then NZC >= IU-IL+1.
216*> If NZC = -1, then a workspace query is assumed; the
217*> routine calculates the number of columns of the array Z that
218*> are needed to hold the eigenvectors.
219*> This value is returned as the first entry of the Z array, and
220*> no error message related to NZC is issued by XERBLA.
221*> \endverbatim
222*>
223*> \param[out] ISUPPZ
224*> \verbatim
225*> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
226*> The support of the eigenvectors in Z, i.e., the indices
227*> indicating the nonzero elements in Z. The i-th computed eigenvector
228*> is nonzero only in elements ISUPPZ( 2*i-1 ) through
229*> ISUPPZ( 2*i ). This is relevant in the case when the matrix
230*> is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
231*> \endverbatim
232*>
233*> \param[in,out] TRYRAC
234*> \verbatim
235*> TRYRAC is LOGICAL
236*> If TRYRAC = .TRUE., indicates that the code should check whether
237*> the tridiagonal matrix defines its eigenvalues to high relative
238*> accuracy. If so, the code uses relative-accuracy preserving
239*> algorithms that might be (a bit) slower depending on the matrix.
240*> If the matrix does not define its eigenvalues to high relative
241*> accuracy, the code can uses possibly faster algorithms.
242*> If TRYRAC = .FALSE., the code is not required to guarantee
243*> relatively accurate eigenvalues and can use the fastest possible
244*> techniques.
245*> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
246*> does not define its eigenvalues to high relative accuracy.
247*> \endverbatim
248*>
249*> \param[out] WORK
250*> \verbatim
251*> WORK is REAL array, dimension (LWORK)
252*> On exit, if INFO = 0, WORK(1) returns the optimal
253*> (and minimal) LWORK.
254*> \endverbatim
255*>
256*> \param[in] LWORK
257*> \verbatim
258*> LWORK is INTEGER
259*> The dimension of the array WORK. LWORK >= max(1,18*N)
260*> if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
261*> If LWORK = -1, then a workspace query is assumed; the routine
262*> only calculates the optimal size of the WORK array, returns
263*> this value as the first entry of the WORK array, and no error
264*> message related to LWORK is issued by XERBLA.
265*> \endverbatim
266*>
267*> \param[out] IWORK
268*> \verbatim
269*> IWORK is INTEGER array, dimension (LIWORK)
270*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
271*> \endverbatim
272*>
273*> \param[in] LIWORK
274*> \verbatim
275*> LIWORK is INTEGER
276*> The dimension of the array IWORK. LIWORK >= max(1,10*N)
277*> if the eigenvectors are desired, and LIWORK >= max(1,8*N)
278*> if only the eigenvalues are to be computed.
279*> If LIWORK = -1, then a workspace query is assumed; the
280*> routine only calculates the optimal size of the IWORK array,
281*> returns this value as the first entry of the IWORK array, and
282*> no error message related to LIWORK is issued by XERBLA.
283*> \endverbatim
284*>
285*> \param[out] INFO
286*> \verbatim
287*> INFO is INTEGER
288*> On exit, INFO
289*> = 0: successful exit
290*> < 0: if INFO = -i, the i-th argument had an illegal value
291*> > 0: if INFO = 1X, internal error in SLARRE,
292*> if INFO = 2X, internal error in SLARRV.
293*> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
294*> the nonzero error code returned by SLARRE or
295*> SLARRV, respectively.
296*> \endverbatim
297*
298* Authors:
299* ========
300*
301*> \author Univ. of Tennessee
302*> \author Univ. of California Berkeley
303*> \author Univ. of Colorado Denver
304*> \author NAG Ltd.
305*
306*> \ingroup stemr
307*
308*> \par Contributors:
309* ==================
310*>
311*> Beresford Parlett, University of California, Berkeley, USA \n
312*> Jim Demmel, University of California, Berkeley, USA \n
313*> Inderjit Dhillon, University of Texas, Austin, USA \n
314*> Osni Marques, LBNL/NERSC, USA \n
315*> Christof Voemel, University of California, Berkeley, USA \n
316*> Aravindh Krishnamoorthy, FAU, Erlangen, Germany \n
317*
318* =====================================================================
319 SUBROUTINE sstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
320 $ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
321 $ IWORK, LIWORK, INFO )
322*
323* -- LAPACK computational routine --
324* -- LAPACK is a software package provided by Univ. of Tennessee, --
325* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
326*
327* .. Scalar Arguments ..
328 CHARACTER JOBZ, RANGE
329 LOGICAL TRYRAC
330 INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
331 REAL VL, VU
332* ..
333* .. Array Arguments ..
334 INTEGER ISUPPZ( * ), IWORK( * )
335 REAL D( * ), E( * ), W( * ), WORK( * )
336 REAL Z( LDZ, * )
337* ..
338*
339* =====================================================================
340*
341* .. Parameters ..
342 REAL ZERO, ONE, FOUR, MINRGP
343 PARAMETER ( ZERO = 0.0e0, one = 1.0e0,
344 $ four = 4.0e0,
345 $ minrgp = 3.0e-3 )
346* ..
347* .. Local Scalars ..
348 LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
349 $ LAESWAP
350 INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
351 $ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
352 $ inde2, inderr, indgp, indgrs, indwrk, itmp,
353 $ itmp2, j, jblk, jj, liwmin, lwmin, nsplit,
354 $ nzcmin, offset, wbegin, wend
355 REAL BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
356 $ RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
357 $ thresh, tmp, tnrm, wl, wu
358* ..
359* ..
360* .. External Functions ..
361 LOGICAL LSAME
362 REAL SLAMCH, SLANST, SROUNDUP_LWORK
363 EXTERNAL lsame, slamch, slanst, sroundup_lwork
364* ..
365* .. External Subroutines ..
366 EXTERNAL scopy, slae2, slaev2, slarrc, slarre, slarrj,
368* ..
369* .. Intrinsic Functions ..
370 INTRINSIC max, min, sqrt
371* ..
372* .. Executable Statements ..
373*
374* Test the input parameters.
375*
376 wantz = lsame( jobz, 'V' )
377 alleig = lsame( range, 'A' )
378 valeig = lsame( range, 'V' )
379 indeig = lsame( range, 'I' )
380*
381 lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
382 zquery = ( nzc.EQ.-1 )
383 laeswap = .false.
384
385* SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
386* In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
387* Furthermore, SLARRV needs WORK of size 12*N, IWORK of size 7*N.
388 IF( wantz ) THEN
389 lwmin = 18*n
390 liwmin = 10*n
391 ELSE
392* need less workspace if only the eigenvalues are wanted
393 lwmin = 12*n
394 liwmin = 8*n
395 ENDIF
396
397 wl = zero
398 wu = zero
399 iil = 0
400 iiu = 0
401 nsplit = 0
402
403 IF( valeig ) THEN
404* We do not reference VL, VU in the cases RANGE = 'I','A'
405* The interval (WL, WU] contains all the wanted eigenvalues.
406* It is either given by the user or computed in SLARRE.
407 wl = vl
408 wu = vu
409 ELSEIF( indeig ) THEN
410* We do not reference IL, IU in the cases RANGE = 'V','A'
411 iil = il
412 iiu = iu
413 ENDIF
414*
415 info = 0
416 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
417 info = -1
418 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
419 info = -2
420 ELSE IF( n.LT.0 ) THEN
421 info = -3
422 ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
423 info = -7
424 ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
425 info = -8
426 ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
427 info = -9
428 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
429 info = -13
430 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
431 info = -17
432 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
433 info = -19
434 END IF
435*
436* Get machine constants.
437*
438 safmin = slamch( 'Safe minimum' )
439 eps = slamch( 'Precision' )
440 smlnum = safmin / eps
441 bignum = one / smlnum
442 rmin = sqrt( smlnum )
443 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
444*
445 IF( info.EQ.0 ) THEN
446 work( 1 ) = sroundup_lwork(lwmin)
447 iwork( 1 ) = liwmin
448*
449 IF( wantz .AND. alleig ) THEN
450 nzcmin = n
451 ELSE IF( wantz .AND. valeig ) THEN
452 CALL slarrc( 'T', n, vl, vu, d, e, safmin,
453 $ nzcmin, itmp, itmp2, info )
454 ELSE IF( wantz .AND. indeig ) THEN
455 nzcmin = iiu-iil+1
456 ELSE
457* WANTZ .EQ. FALSE.
458 nzcmin = 0
459 ENDIF
460 IF( zquery .AND. info.EQ.0 ) THEN
461 z( 1,1 ) = nzcmin
462 ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
463 info = -14
464 END IF
465 END IF
466
467 IF( info.NE.0 ) THEN
468*
469 CALL xerbla( 'SSTEMR', -info )
470*
471 RETURN
472 ELSE IF( lquery .OR. zquery ) THEN
473 RETURN
474 END IF
475*
476* Handle N = 0, 1, and 2 cases immediately
477*
478 m = 0
479 IF( n.EQ.0 )
480 $ RETURN
481*
482 IF( n.EQ.1 ) THEN
483 IF( alleig .OR. indeig ) THEN
484 m = 1
485 w( 1 ) = d( 1 )
486 ELSE
487 IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
488 m = 1
489 w( 1 ) = d( 1 )
490 END IF
491 END IF
492 IF( wantz.AND.(.NOT.zquery) ) THEN
493 z( 1, 1 ) = one
494 isuppz(1) = 1
495 isuppz(2) = 1
496 END IF
497 RETURN
498 END IF
499*
500 IF( n.EQ.2 ) THEN
501 IF( .NOT.wantz ) THEN
502 CALL slae2( d(1), e(1), d(2), r1, r2 )
503 ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
504 CALL slaev2( d(1), e(1), d(2), r1, r2, cs, sn )
505 END IF
506* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
507* the following code requires R1 >= R2. Hence, we correct
508* the order of R1, R2, CS, SN if R1 < R2 before further processing.
509 IF( r1.LT.r2 ) THEN
510 e(2) = r1
511 r1 = r2
512 r2 = e(2)
513 laeswap = .true.
514 ENDIF
515 IF( alleig.OR.
516 $ (valeig.AND.(r2.GT.wl).AND.
517 $ (r2.LE.wu)).OR.
518 $ (indeig.AND.(iil.EQ.1)) ) THEN
519 m = m+1
520 w( m ) = r2
521 IF( wantz.AND.(.NOT.zquery) ) THEN
522 IF( laeswap ) THEN
523 z( 1, m ) = cs
524 z( 2, m ) = sn
525 ELSE
526 z( 1, m ) = -sn
527 z( 2, m ) = cs
528 ENDIF
529* Note: At most one of SN and CS can be zero.
530 IF (sn.NE.zero) THEN
531 IF (cs.NE.zero) THEN
532 isuppz(2*m-1) = 1
533 isuppz(2*m) = 2
534 ELSE
535 isuppz(2*m-1) = 1
536 isuppz(2*m) = 1
537 END IF
538 ELSE
539 isuppz(2*m-1) = 2
540 isuppz(2*m) = 2
541 END IF
542 ENDIF
543 ENDIF
544 IF( alleig.OR.
545 $ (valeig.AND.(r1.GT.wl).AND.
546 $ (r1.LE.wu)).OR.
547 $ (indeig.AND.(iiu.EQ.2)) ) THEN
548 m = m+1
549 w( m ) = r1
550 IF( wantz.AND.(.NOT.zquery) ) THEN
551 IF( laeswap ) THEN
552 z( 1, m ) = -sn
553 z( 2, m ) = cs
554 ELSE
555 z( 1, m ) = cs
556 z( 2, m ) = sn
557 ENDIF
558* Note: At most one of SN and CS can be zero.
559 IF (sn.NE.zero) THEN
560 IF (cs.NE.zero) THEN
561 isuppz(2*m-1) = 1
562 isuppz(2*m) = 2
563 ELSE
564 isuppz(2*m-1) = 1
565 isuppz(2*m) = 1
566 END IF
567 ELSE
568 isuppz(2*m-1) = 2
569 isuppz(2*m) = 2
570 END IF
571 ENDIF
572 ENDIF
573 ELSE
574
575* Continue with general N
576
577 indgrs = 1
578 inderr = 2*n + 1
579 indgp = 3*n + 1
580 indd = 4*n + 1
581 inde2 = 5*n + 1
582 indwrk = 6*n + 1
583*
584 iinspl = 1
585 iindbl = n + 1
586 iindw = 2*n + 1
587 iindwk = 3*n + 1
588*
589* Scale matrix to allowable range, if necessary.
590* The allowable range is related to the PIVMIN parameter; see the
591* comments in SLARRD. The preference for scaling small values
592* up is heuristic; we expect users' matrices not to be close to the
593* RMAX threshold.
594*
595 scale = one
596 tnrm = slanst( 'M', n, d, e )
597 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
598 scale = rmin / tnrm
599 ELSE IF( tnrm.GT.rmax ) THEN
600 scale = rmax / tnrm
601 END IF
602 IF( scale.NE.one ) THEN
603 CALL sscal( n, scale, d, 1 )
604 CALL sscal( n-1, scale, e, 1 )
605 tnrm = tnrm*scale
606 IF( valeig ) THEN
607* If eigenvalues in interval have to be found,
608* scale (WL, WU] accordingly
609 wl = wl*scale
610 wu = wu*scale
611 ENDIF
612 END IF
613*
614* Compute the desired eigenvalues of the tridiagonal after splitting
615* into smaller subblocks if the corresponding off-diagonal elements
616* are small
617* THRESH is the splitting parameter for SLARRE
618* A negative THRESH forces the old splitting criterion based on the
619* size of the off-diagonal. A positive THRESH switches to splitting
620* which preserves relative accuracy.
621*
622 IF( tryrac ) THEN
623* Test whether the matrix warrants the more expensive relative approach.
624 CALL slarrr( n, d, e, iinfo )
625 ELSE
626* The user does not care about relative accurately eigenvalues
627 iinfo = -1
628 ENDIF
629* Set the splitting criterion
630 IF (iinfo.EQ.0) THEN
631 thresh = eps
632 ELSE
633 thresh = -eps
634* relative accuracy is desired but T does not guarantee it
635 tryrac = .false.
636 ENDIF
637*
638 IF( tryrac ) THEN
639* Copy original diagonal, needed to guarantee relative accuracy
640 CALL scopy(n,d,1,work(indd),1)
641 ENDIF
642* Store the squares of the offdiagonal values of T
643 DO 5 j = 1, n-1
644 work( inde2+j-1 ) = e(j)**2
645 5 CONTINUE
646
647* Set the tolerance parameters for bisection
648 IF( .NOT.wantz ) THEN
649* SLARRE computes the eigenvalues to full precision.
650 rtol1 = four * eps
651 rtol2 = four * eps
652 ELSE
653* SLARRE computes the eigenvalues to less than full precision.
654* SLARRV will refine the eigenvalue approximations, and we can
655* need less accurate initial bisection in SLARRE.
656* Note: these settings do only affect the subset case and SLARRE
657 rtol1 = max( sqrt(eps)*5.0e-2, four * eps )
658 rtol2 = max( sqrt(eps)*5.0e-3, four * eps )
659 ENDIF
660 CALL slarre( range, n, wl, wu, iil, iiu, d, e,
661 $ work(inde2), rtol1, rtol2, thresh, nsplit,
662 $ iwork( iinspl ), m, w, work( inderr ),
663 $ work( indgp ), iwork( iindbl ),
664 $ iwork( iindw ), work( indgrs ), pivmin,
665 $ work( indwrk ), iwork( iindwk ), iinfo )
666 IF( iinfo.NE.0 ) THEN
667 info = 10 + abs( iinfo )
668 RETURN
669 END IF
670* Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
671* part of the spectrum. All desired eigenvalues are contained in
672* (WL,WU]
673
674
675 IF( wantz ) THEN
676*
677* Compute the desired eigenvectors corresponding to the computed
678* eigenvalues
679*
680 CALL slarrv( n, wl, wu, d, e,
681 $ pivmin, iwork( iinspl ), m,
682 $ 1, m, minrgp, rtol1, rtol2,
683 $ w, work( inderr ), work( indgp ), iwork( iindbl ),
684 $ iwork( iindw ), work( indgrs ), z, ldz,
685 $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
686 IF( iinfo.NE.0 ) THEN
687 info = 20 + abs( iinfo )
688 RETURN
689 END IF
690 ELSE
691* SLARRE computes eigenvalues of the (shifted) root representation
692* SLARRV returns the eigenvalues of the unshifted matrix.
693* However, if the eigenvectors are not desired by the user, we need
694* to apply the corresponding shifts from SLARRE to obtain the
695* eigenvalues of the original matrix.
696 DO 20 j = 1, m
697 itmp = iwork( iindbl+j-1 )
698 w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
699 20 CONTINUE
700 END IF
701*
702
703 IF ( tryrac ) THEN
704* Refine computed eigenvalues so that they are relatively accurate
705* with respect to the original matrix T.
706 ibegin = 1
707 wbegin = 1
708 DO 39 jblk = 1, iwork( iindbl+m-1 )
709 iend = iwork( iinspl+jblk-1 )
710 in = iend - ibegin + 1
711 wend = wbegin - 1
712* check if any eigenvalues have to be refined in this block
713 36 CONTINUE
714 IF( wend.LT.m ) THEN
715 IF( iwork( iindbl+wend ).EQ.jblk ) THEN
716 wend = wend + 1
717 GO TO 36
718 END IF
719 END IF
720 IF( wend.LT.wbegin ) THEN
721 ibegin = iend + 1
722 GO TO 39
723 END IF
724
725 offset = iwork(iindw+wbegin-1)-1
726 ifirst = iwork(iindw+wbegin-1)
727 ilast = iwork(iindw+wend-1)
728 rtol2 = four * eps
729 CALL slarrj( in,
730 $ work(indd+ibegin-1), work(inde2+ibegin-1),
731 $ ifirst, ilast, rtol2, offset, w(wbegin),
732 $ work( inderr+wbegin-1 ),
733 $ work( indwrk ), iwork( iindwk ), pivmin,
734 $ tnrm, iinfo )
735 ibegin = iend + 1
736 wbegin = wend + 1
737 39 CONTINUE
738 ENDIF
739*
740* If matrix was scaled, then rescale eigenvalues appropriately.
741*
742 IF( scale.NE.one ) THEN
743 CALL sscal( m, one / scale, w, 1 )
744 END IF
745 END IF
746*
747* If eigenvalues are not in increasing order, then sort them,
748* possibly along with eigenvectors.
749*
750 IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
751 IF( .NOT. wantz ) THEN
752 CALL slasrt( 'I', m, w, iinfo )
753 IF( iinfo.NE.0 ) THEN
754 info = 3
755 RETURN
756 END IF
757 ELSE
758 DO 60 j = 1, m - 1
759 i = 0
760 tmp = w( j )
761 DO 50 jj = j + 1, m
762 IF( w( jj ).LT.tmp ) THEN
763 i = jj
764 tmp = w( jj )
765 END IF
766 50 CONTINUE
767 IF( i.NE.0 ) THEN
768 w( i ) = w( j )
769 w( j ) = tmp
770 IF( wantz ) THEN
771 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
772 itmp = isuppz( 2*i-1 )
773 isuppz( 2*i-1 ) = isuppz( 2*j-1 )
774 isuppz( 2*j-1 ) = itmp
775 itmp = isuppz( 2*i )
776 isuppz( 2*i ) = isuppz( 2*j )
777 isuppz( 2*j ) = itmp
778 END IF
779 END IF
780 60 CONTINUE
781 END IF
782 ENDIF
783*
784*
785 work( 1 ) = sroundup_lwork(lwmin)
786 iwork( 1 ) = liwmin
787 RETURN
788*
789* End of SSTEMR
790*
791 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine slae2(a, b, c, rt1, rt2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition slae2.f:102
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:120
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:137
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:305
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:168
subroutine slarrr(n, d, e, info)
SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computa...
Definition slarrr.f:94
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:292
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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:322
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82