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