LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dstemr.f
Go to the documentation of this file.
1*> \brief \b DSTEMR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSTEMR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstemr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstemr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstemr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSTEMR( 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* DOUBLE PRECISION VL, VU
30* ..
31* .. Array Arguments ..
32* INTEGER ISUPPZ( * ), IWORK( * )
33* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
34* DOUBLE PRECISION Z( LDZ, * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> DSTEMR 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.DSTEMR 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLARRE,
292*> if INFO = 2X, internal error in DLARRV.
293*> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
294*> the nonzero error code returned by DLARRE or
295*> DLARRV, 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 dstemr( 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 DOUBLE PRECISION VL, VU
332* ..
333* .. Array Arguments ..
334 INTEGER ISUPPZ( * ), IWORK( * )
335 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
336 DOUBLE PRECISION Z( LDZ, * )
337* ..
338*
339* =====================================================================
340*
341* .. Parameters ..
342 DOUBLE PRECISION ZERO, ONE, FOUR, MINRGP
343 PARAMETER ( ZERO = 0.0d0, one = 1.0d0,
344 $ four = 4.0d0,
345 $ minrgp = 1.0d-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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH, DLANST
363 EXTERNAL lsame, dlamch, dlanst
364* ..
365* .. External Subroutines ..
366 EXTERNAL dcopy, dlae2, dlaev2, dlarrc, dlarre, dlarrj,
368* ..
369* .. Intrinsic Functions ..
370 INTRINSIC max, min, sqrt
371
372
373* ..
374* .. Executable Statements ..
375*
376* Test the input parameters.
377*
378 wantz = lsame( jobz, 'V' )
379 alleig = lsame( range, 'A' )
380 valeig = lsame( range, 'V' )
381 indeig = lsame( range, 'I' )
382*
383 lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
384 zquery = ( nzc.EQ.-1 )
385 laeswap = .false.
386
387* DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
388* In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
389* Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N.
390 IF( wantz ) THEN
391 lwmin = 18*n
392 liwmin = 10*n
393 ELSE
394* need less workspace if only the eigenvalues are wanted
395 lwmin = 12*n
396 liwmin = 8*n
397 ENDIF
398
399 wl = zero
400 wu = zero
401 iil = 0
402 iiu = 0
403 nsplit = 0
404
405 IF( valeig ) THEN
406* We do not reference VL, VU in the cases RANGE = 'I','A'
407* The interval (WL, WU] contains all the wanted eigenvalues.
408* It is either given by the user or computed in DLARRE.
409 wl = vl
410 wu = vu
411 ELSEIF( indeig ) THEN
412* We do not reference IL, IU in the cases RANGE = 'V','A'
413 iil = il
414 iiu = iu
415 ENDIF
416*
417 info = 0
418 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
419 info = -1
420 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
421 info = -2
422 ELSE IF( n.LT.0 ) THEN
423 info = -3
424 ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
425 info = -7
426 ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
427 info = -8
428 ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
429 info = -9
430 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
431 info = -13
432 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
433 info = -17
434 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
435 info = -19
436 END IF
437*
438* Get machine constants.
439*
440 safmin = dlamch( 'Safe minimum' )
441 eps = dlamch( 'Precision' )
442 smlnum = safmin / eps
443 bignum = one / smlnum
444 rmin = sqrt( smlnum )
445 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
446*
447 IF( info.EQ.0 ) THEN
448 work( 1 ) = lwmin
449 iwork( 1 ) = liwmin
450*
451 IF( wantz .AND. alleig ) THEN
452 nzcmin = n
453 ELSE IF( wantz .AND. valeig ) THEN
454 CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
455 $ nzcmin, itmp, itmp2, info )
456 ELSE IF( wantz .AND. indeig ) THEN
457 nzcmin = iiu-iil+1
458 ELSE
459* WANTZ .EQ. FALSE.
460 nzcmin = 0
461 ENDIF
462 IF( zquery .AND. info.EQ.0 ) THEN
463 z( 1,1 ) = nzcmin
464 ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
465 info = -14
466 END IF
467 END IF
468
469 IF( info.NE.0 ) THEN
470*
471 CALL xerbla( 'DSTEMR', -info )
472*
473 RETURN
474 ELSE IF( lquery .OR. zquery ) THEN
475 RETURN
476 END IF
477*
478* Handle N = 0, 1, and 2 cases immediately
479*
480 m = 0
481 IF( n.EQ.0 )
482 $ RETURN
483*
484 IF( n.EQ.1 ) THEN
485 IF( alleig .OR. indeig ) THEN
486 m = 1
487 w( 1 ) = d( 1 )
488 ELSE
489 IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
490 m = 1
491 w( 1 ) = d( 1 )
492 END IF
493 END IF
494 IF( wantz.AND.(.NOT.zquery) ) THEN
495 z( 1, 1 ) = one
496 isuppz(1) = 1
497 isuppz(2) = 1
498 END IF
499 RETURN
500 END IF
501*
502 IF( n.EQ.2 ) THEN
503 IF( .NOT.wantz ) THEN
504 CALL dlae2( d(1), e(1), d(2), r1, r2 )
505 ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
506 CALL dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
507 END IF
508* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
509* the following code requires R1 >= R2. Hence, we correct
510* the order of R1, R2, CS, SN if R1 < R2 before further processing.
511 IF( r1.LT.r2 ) THEN
512 e(2) = r1
513 r1 = r2
514 r2 = e(2)
515 laeswap = .true.
516 ENDIF
517 IF( alleig.OR.
518 $ (valeig.AND.(r2.GT.wl).AND.
519 $ (r2.LE.wu)).OR.
520 $ (indeig.AND.(iil.EQ.1)) ) THEN
521 m = m+1
522 w( m ) = r2
523 IF( wantz.AND.(.NOT.zquery) ) THEN
524 IF( laeswap ) THEN
525 z( 1, m ) = cs
526 z( 2, m ) = sn
527 ELSE
528 z( 1, m ) = -sn
529 z( 2, m ) = cs
530 ENDIF
531* Note: At most one of SN and CS can be zero.
532 IF (sn.NE.zero) THEN
533 IF (cs.NE.zero) THEN
534 isuppz(2*m-1) = 1
535 isuppz(2*m) = 2
536 ELSE
537 isuppz(2*m-1) = 1
538 isuppz(2*m) = 1
539 END IF
540 ELSE
541 isuppz(2*m-1) = 2
542 isuppz(2*m) = 2
543 END IF
544 ENDIF
545 ENDIF
546 IF( alleig.OR.
547 $ (valeig.AND.(r1.GT.wl).AND.
548 $ (r1.LE.wu)).OR.
549 $ (indeig.AND.(iiu.EQ.2)) ) THEN
550 m = m+1
551 w( m ) = r1
552 IF( wantz.AND.(.NOT.zquery) ) THEN
553 IF( laeswap ) THEN
554 z( 1, m ) = -sn
555 z( 2, m ) = cs
556 ELSE
557 z( 1, m ) = cs
558 z( 2, m ) = sn
559 ENDIF
560* Note: At most one of SN and CS can be zero.
561 IF (sn.NE.zero) THEN
562 IF (cs.NE.zero) THEN
563 isuppz(2*m-1) = 1
564 isuppz(2*m) = 2
565 ELSE
566 isuppz(2*m-1) = 1
567 isuppz(2*m) = 1
568 END IF
569 ELSE
570 isuppz(2*m-1) = 2
571 isuppz(2*m) = 2
572 END IF
573 ENDIF
574 ENDIF
575
576 ELSE
577
578* Continue with general N
579
580 indgrs = 1
581 inderr = 2*n + 1
582 indgp = 3*n + 1
583 indd = 4*n + 1
584 inde2 = 5*n + 1
585 indwrk = 6*n + 1
586*
587 iinspl = 1
588 iindbl = n + 1
589 iindw = 2*n + 1
590 iindwk = 3*n + 1
591*
592* Scale matrix to allowable range, if necessary.
593* The allowable range is related to the PIVMIN parameter; see the
594* comments in DLARRD. The preference for scaling small values
595* up is heuristic; we expect users' matrices not to be close to the
596* RMAX threshold.
597*
598 scale = one
599 tnrm = dlanst( 'M', n, d, e )
600 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
601 scale = rmin / tnrm
602 ELSE IF( tnrm.GT.rmax ) THEN
603 scale = rmax / tnrm
604 END IF
605 IF( scale.NE.one ) THEN
606 CALL dscal( n, scale, d, 1 )
607 CALL dscal( n-1, scale, e, 1 )
608 tnrm = tnrm*scale
609 IF( valeig ) THEN
610* If eigenvalues in interval have to be found,
611* scale (WL, WU] accordingly
612 wl = wl*scale
613 wu = wu*scale
614 ENDIF
615 END IF
616*
617* Compute the desired eigenvalues of the tridiagonal after splitting
618* into smaller subblocks if the corresponding off-diagonal elements
619* are small
620* THRESH is the splitting parameter for DLARRE
621* A negative THRESH forces the old splitting criterion based on the
622* size of the off-diagonal. A positive THRESH switches to splitting
623* which preserves relative accuracy.
624*
625 IF( tryrac ) THEN
626* Test whether the matrix warrants the more expensive relative approach.
627 CALL dlarrr( n, d, e, iinfo )
628 ELSE
629* The user does not care about relative accurately eigenvalues
630 iinfo = -1
631 ENDIF
632* Set the splitting criterion
633 IF (iinfo.EQ.0) THEN
634 thresh = eps
635 ELSE
636 thresh = -eps
637* relative accuracy is desired but T does not guarantee it
638 tryrac = .false.
639 ENDIF
640*
641 IF( tryrac ) THEN
642* Copy original diagonal, needed to guarantee relative accuracy
643 CALL dcopy(n,d,1,work(indd),1)
644 ENDIF
645* Store the squares of the offdiagonal values of T
646 DO 5 j = 1, n-1
647 work( inde2+j-1 ) = e(j)**2
648 5 CONTINUE
649
650* Set the tolerance parameters for bisection
651 IF( .NOT.wantz ) THEN
652* DLARRE computes the eigenvalues to full precision.
653 rtol1 = four * eps
654 rtol2 = four * eps
655 ELSE
656* DLARRE computes the eigenvalues to less than full precision.
657* DLARRV will refine the eigenvalue approximations, and we can
658* need less accurate initial bisection in DLARRE.
659* Note: these settings do only affect the subset case and DLARRE
660 rtol1 = sqrt(eps)
661 rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
662 ENDIF
663 CALL dlarre( range, n, wl, wu, iil, iiu, d, e,
664 $ work(inde2), rtol1, rtol2, thresh, nsplit,
665 $ iwork( iinspl ), m, w, work( inderr ),
666 $ work( indgp ), iwork( iindbl ),
667 $ iwork( iindw ), work( indgrs ), pivmin,
668 $ work( indwrk ), iwork( iindwk ), iinfo )
669 IF( iinfo.NE.0 ) THEN
670 info = 10 + abs( iinfo )
671 RETURN
672 END IF
673* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
674* part of the spectrum. All desired eigenvalues are contained in
675* (WL,WU]
676
677
678 IF( wantz ) THEN
679*
680* Compute the desired eigenvectors corresponding to the computed
681* eigenvalues
682*
683 CALL dlarrv( n, wl, wu, d, e,
684 $ pivmin, iwork( iinspl ), m,
685 $ 1, m, minrgp, rtol1, rtol2,
686 $ w, work( inderr ), work( indgp ), iwork( iindbl ),
687 $ iwork( iindw ), work( indgrs ), z, ldz,
688 $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
689 IF( iinfo.NE.0 ) THEN
690 info = 20 + abs( iinfo )
691 RETURN
692 END IF
693 ELSE
694* DLARRE computes eigenvalues of the (shifted) root representation
695* DLARRV returns the eigenvalues of the unshifted matrix.
696* However, if the eigenvectors are not desired by the user, we need
697* to apply the corresponding shifts from DLARRE to obtain the
698* eigenvalues of the original matrix.
699 DO 20 j = 1, m
700 itmp = iwork( iindbl+j-1 )
701 w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
702 20 CONTINUE
703 END IF
704*
705
706 IF ( tryrac ) THEN
707* Refine computed eigenvalues so that they are relatively accurate
708* with respect to the original matrix T.
709 ibegin = 1
710 wbegin = 1
711 DO 39 jblk = 1, iwork( iindbl+m-1 )
712 iend = iwork( iinspl+jblk-1 )
713 in = iend - ibegin + 1
714 wend = wbegin - 1
715* check if any eigenvalues have to be refined in this block
716 36 CONTINUE
717 IF( wend.LT.m ) THEN
718 IF( iwork( iindbl+wend ).EQ.jblk ) THEN
719 wend = wend + 1
720 GO TO 36
721 END IF
722 END IF
723 IF( wend.LT.wbegin ) THEN
724 ibegin = iend + 1
725 GO TO 39
726 END IF
727
728 offset = iwork(iindw+wbegin-1)-1
729 ifirst = iwork(iindw+wbegin-1)
730 ilast = iwork(iindw+wend-1)
731 rtol2 = four * eps
732 CALL dlarrj( in,
733 $ work(indd+ibegin-1), work(inde2+ibegin-1),
734 $ ifirst, ilast, rtol2, offset, w(wbegin),
735 $ work( inderr+wbegin-1 ),
736 $ work( indwrk ), iwork( iindwk ), pivmin,
737 $ tnrm, iinfo )
738 ibegin = iend + 1
739 wbegin = wend + 1
740 39 CONTINUE
741 ENDIF
742*
743* If matrix was scaled, then rescale eigenvalues appropriately.
744*
745 IF( scale.NE.one ) THEN
746 CALL dscal( m, one / scale, w, 1 )
747 END IF
748
749 END IF
750
751*
752* If eigenvalues are not in increasing order, then sort them,
753* possibly along with eigenvectors.
754*
755 IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
756 IF( .NOT. wantz ) THEN
757 CALL dlasrt( 'I', m, w, iinfo )
758 IF( iinfo.NE.0 ) THEN
759 info = 3
760 RETURN
761 END IF
762 ELSE
763 DO 60 j = 1, m - 1
764 i = 0
765 tmp = w( j )
766 DO 50 jj = j + 1, m
767 IF( w( jj ).LT.tmp ) THEN
768 i = jj
769 tmp = w( jj )
770 END IF
771 50 CONTINUE
772 IF( i.NE.0 ) THEN
773 w( i ) = w( j )
774 w( j ) = tmp
775 IF( wantz ) THEN
776 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
777 itmp = isuppz( 2*i-1 )
778 isuppz( 2*i-1 ) = isuppz( 2*j-1 )
779 isuppz( 2*j-1 ) = itmp
780 itmp = isuppz( 2*i )
781 isuppz( 2*i ) = isuppz( 2*j )
782 isuppz( 2*j ) = itmp
783 END IF
784 END IF
785 60 CONTINUE
786 END IF
787 ENDIF
788*
789*
790 work( 1 ) = lwmin
791 iwork( 1 ) = liwmin
792 RETURN
793*
794* End of DSTEMR
795*
796 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
subroutine dlaev2(a, b, c, rt1, rt2, cs1, sn1)
DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition dlaev2.f:120
subroutine dlarrc(jobt, n, vl, vu, d, e, pivmin, eigcnt, lcnt, rcnt, info)
DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
Definition dlarrc.f:137
subroutine dlarre(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)
DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
Definition dlarre.f:305
subroutine dlarrj(n, d, e2, ifirst, ilast, rtol, offset, w, werr, work, iwork, pivmin, spdiam, info)
DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
Definition dlarrj.f:168
subroutine dlarrr(n, d, e, info)
DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computa...
Definition dlarrr.f:94
subroutine dlarrv(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)
DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition dlarrv.f:292
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
DSTEMR
Definition dstemr.f:322
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82