LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dstevr.f
Go to the documentation of this file.
1*> \brief <b> DSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstevr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstevr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstevr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
22* M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
23* LIWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER JOBZ, RANGE
27* INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
28* DOUBLE PRECISION ABSTOL, VL, VU
29* ..
30* .. Array Arguments ..
31* INTEGER ISUPPZ( * ), IWORK( * )
32* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> DSTEVR computes selected eigenvalues and, optionally, eigenvectors
42*> of a real symmetric tridiagonal matrix T. Eigenvalues and
43*> eigenvectors can be selected by specifying either a range of values
44*> or a range of indices for the desired eigenvalues.
45*>
46*> Whenever possible, DSTEVR calls DSTEMR to compute the
47*> eigenspectrum using Relatively Robust Representations. DSTEMR
48*> computes eigenvalues by the dqds algorithm, while orthogonal
49*> eigenvectors are computed from various "good" L D L^T representations
50*> (also known as Relatively Robust Representations). Gram-Schmidt
51*> orthogonalization is avoided as far as possible. More specifically,
52*> the various steps of the algorithm are as follows. For the i-th
53*> unreduced block of T,
54*> (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
55*> is a relatively robust representation,
56*> (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
57*> relative accuracy by the dqds algorithm,
58*> (c) If there is a cluster of close eigenvalues, "choose" sigma_i
59*> close to the cluster, and go to step (a),
60*> (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
61*> compute the corresponding eigenvector by forming a
62*> rank-revealing twisted factorization.
63*> The desired accuracy of the output can be specified by the input
64*> parameter ABSTOL.
65*>
66*> For more details, see "A new O(n^2) algorithm for the symmetric
67*> tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
68*> Computer Science Division Technical Report No. UCB//CSD-97-971,
69*> UC Berkeley, May 1997.
70*>
71*>
72*> Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested
73*> on machines which conform to the ieee-754 floating point standard.
74*> DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and
75*> when partial spectrum requests are made.
76*>
77*> Normal execution of DSTEMR may create NaNs and infinities and
78*> hence may abort due to a floating point exception in environments
79*> which do not handle NaNs and infinities in the ieee standard default
80*> manner.
81*> \endverbatim
82*
83* Arguments:
84* ==========
85*
86*> \param[in] JOBZ
87*> \verbatim
88*> JOBZ is CHARACTER*1
89*> = 'N': Compute eigenvalues only;
90*> = 'V': Compute eigenvalues and eigenvectors.
91*> \endverbatim
92*>
93*> \param[in] RANGE
94*> \verbatim
95*> RANGE is CHARACTER*1
96*> = 'A': all eigenvalues will be found.
97*> = 'V': all eigenvalues in the half-open interval (VL,VU]
98*> will be found.
99*> = 'I': the IL-th through IU-th eigenvalues will be found.
100*> For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
101*> DSTEIN are called
102*> \endverbatim
103*>
104*> \param[in] N
105*> \verbatim
106*> N is INTEGER
107*> The order of the matrix. N >= 0.
108*> \endverbatim
109*>
110*> \param[in,out] D
111*> \verbatim
112*> D is DOUBLE PRECISION array, dimension (N)
113*> On entry, the n diagonal elements of the tridiagonal matrix
114*> A.
115*> On exit, D may be multiplied by a constant factor chosen
116*> to avoid over/underflow in computing the eigenvalues.
117*> \endverbatim
118*>
119*> \param[in,out] E
120*> \verbatim
121*> E is DOUBLE PRECISION array, dimension (max(1,N-1))
122*> On entry, the (n-1) subdiagonal elements of the tridiagonal
123*> matrix A in elements 1 to N-1 of E.
124*> On exit, E may be multiplied by a constant factor chosen
125*> to avoid over/underflow in computing the eigenvalues.
126*> \endverbatim
127*>
128*> \param[in] VL
129*> \verbatim
130*> VL is DOUBLE PRECISION
131*> If RANGE='V', the lower bound of the interval to
132*> be searched for eigenvalues. VL < VU.
133*> Not referenced if RANGE = 'A' or 'I'.
134*> \endverbatim
135*>
136*> \param[in] VU
137*> \verbatim
138*> VU is DOUBLE PRECISION
139*> If RANGE='V', the upper bound of the interval to
140*> be searched for eigenvalues. VL < VU.
141*> Not referenced if RANGE = 'A' or 'I'.
142*> \endverbatim
143*>
144*> \param[in] IL
145*> \verbatim
146*> IL is INTEGER
147*> If RANGE='I', the index of the
148*> smallest eigenvalue to be returned.
149*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
150*> Not referenced if RANGE = 'A' or 'V'.
151*> \endverbatim
152*>
153*> \param[in] IU
154*> \verbatim
155*> IU is INTEGER
156*> If RANGE='I', the index of the
157*> largest eigenvalue to be returned.
158*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
159*> Not referenced if RANGE = 'A' or 'V'.
160*> \endverbatim
161*>
162*> \param[in] ABSTOL
163*> \verbatim
164*> ABSTOL is DOUBLE PRECISION
165*> The absolute error tolerance for the eigenvalues.
166*> An approximate eigenvalue is accepted as converged
167*> when it is determined to lie in an interval [a,b]
168*> of width less than or equal to
169*>
170*> ABSTOL + EPS * max( |a|,|b| ) ,
171*>
172*> where EPS is the machine precision. If ABSTOL is less than
173*> or equal to zero, then EPS*|T| will be used in its place,
174*> where |T| is the 1-norm of the tridiagonal matrix obtained
175*> by reducing A to tridiagonal form.
176*>
177*> See "Computing Small Singular Values of Bidiagonal Matrices
178*> with Guaranteed High Relative Accuracy," by Demmel and
179*> Kahan, LAPACK Working Note #3.
180*>
181*> If high relative accuracy is important, set ABSTOL to
182*> DLAMCH( 'Safe minimum' ). Doing so will guarantee that
183*> eigenvalues are computed to high relative accuracy when
184*> possible in future releases. The current code does not
185*> make any guarantees about high relative accuracy, but
186*> future releases will. See J. Barlow and J. Demmel,
187*> "Computing Accurate Eigensystems of Scaled Diagonally
188*> Dominant Matrices", LAPACK Working Note #7, for a discussion
189*> of which matrices define their eigenvalues to high relative
190*> accuracy.
191*> \endverbatim
192*>
193*> \param[out] M
194*> \verbatim
195*> M is INTEGER
196*> The total number of eigenvalues found. 0 <= M <= N.
197*> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
198*> \endverbatim
199*>
200*> \param[out] W
201*> \verbatim
202*> W is DOUBLE PRECISION array, dimension (N)
203*> The first M elements contain the selected eigenvalues in
204*> ascending order.
205*> \endverbatim
206*>
207*> \param[out] Z
208*> \verbatim
209*> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
210*> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
211*> contain the orthonormal eigenvectors of the matrix A
212*> corresponding to the selected eigenvalues, with the i-th
213*> column of Z holding the eigenvector associated with W(i).
214*> Note: the user must ensure that at least max(1,M) columns are
215*> supplied in the array Z; if RANGE = 'V', the exact value of M
216*> is not known in advance and an upper bound must be used.
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', LDZ >= max(1,N).
224*> \endverbatim
225*>
226*> \param[out] ISUPPZ
227*> \verbatim
228*> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
229*> The support of the eigenvectors in Z, i.e., the indices
230*> indicating the nonzero elements in Z. The i-th eigenvector
231*> is nonzero only in elements ISUPPZ( 2*i-1 ) through
232*> ISUPPZ( 2*i ).
233*> Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
234*> \endverbatim
235*>
236*> \param[out] WORK
237*> \verbatim
238*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
239*> On exit, if INFO = 0, WORK(1) returns the optimal (and
240*> minimal) LWORK.
241*> \endverbatim
242*>
243*> \param[in] LWORK
244*> \verbatim
245*> LWORK is INTEGER
246*> The dimension of the array WORK. LWORK >= max(1,20*N).
247*>
248*> If LWORK = -1, then a workspace query is assumed; the routine
249*> only calculates the optimal sizes of the WORK and IWORK
250*> arrays, returns these values as the first entries of the WORK
251*> and IWORK arrays, and no error message related to LWORK or
252*> LIWORK is issued by XERBLA.
253*> \endverbatim
254*>
255*> \param[out] IWORK
256*> \verbatim
257*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
258*> On exit, if INFO = 0, IWORK(1) returns the optimal (and
259*> minimal) LIWORK.
260*> \endverbatim
261*>
262*> \param[in] LIWORK
263*> \verbatim
264*> LIWORK is INTEGER
265*> The dimension of the array IWORK. LIWORK >= max(1,10*N).
266*>
267*> If LIWORK = -1, then a workspace query is assumed; the
268*> routine only calculates the optimal sizes of the WORK and
269*> IWORK arrays, returns these values as the first entries of
270*> the WORK and IWORK arrays, and no error message related to
271*> LWORK or LIWORK is issued by XERBLA.
272*> \endverbatim
273*>
274*> \param[out] INFO
275*> \verbatim
276*> INFO is INTEGER
277*> = 0: successful exit
278*> < 0: if INFO = -i, the i-th argument had an illegal value
279*> > 0: Internal error
280*> \endverbatim
281*
282* Authors:
283* ========
284*
285*> \author Univ. of Tennessee
286*> \author Univ. of California Berkeley
287*> \author Univ. of Colorado Denver
288*> \author NAG Ltd.
289*
290*> \ingroup stevr
291*
292*> \par Contributors:
293* ==================
294*>
295*> Inderjit Dhillon, IBM Almaden, USA \n
296*> Osni Marques, LBNL/NERSC, USA \n
297*> Ken Stanley, Computer Science Division, University of
298*> California at Berkeley, USA \n
299*>
300* =====================================================================
301 SUBROUTINE dstevr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
302 \$ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
303 \$ LIWORK, INFO )
304*
305* -- LAPACK driver routine --
306* -- LAPACK is a software package provided by Univ. of Tennessee, --
307* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
308*
309* .. Scalar Arguments ..
310 CHARACTER JOBZ, RANGE
311 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
312 DOUBLE PRECISION ABSTOL, VL, VU
313* ..
314* .. Array Arguments ..
315 INTEGER ISUPPZ( * ), IWORK( * )
316 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
317* ..
318*
319* =====================================================================
320*
321* .. Parameters ..
322 DOUBLE PRECISION ZERO, ONE, TWO
323 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
324* ..
325* .. Local Scalars ..
326 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ,
327 \$ TRYRAC
328 CHARACTER ORDER
329 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP,
330 \$ indiwo, iscale, itmp1, j, jj, liwmin, lwmin,
331 \$ nsplit
332 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
333 \$ TMP1, TNRM, VLL, VUU
334* ..
335* .. External Functions ..
336 LOGICAL LSAME
337 INTEGER ILAENV
338 DOUBLE PRECISION DLAMCH, DLANST
339 EXTERNAL lsame, ilaenv, dlamch, dlanst
340* ..
341* .. External Subroutines ..
342 EXTERNAL dcopy, dscal, dstebz, dstemr, dstein, dsterf,
343 \$ dswap, xerbla
344* ..
345* .. Intrinsic Functions ..
346 INTRINSIC max, min, sqrt
347* ..
348* .. Executable Statements ..
349*
350*
351* Test the input parameters.
352*
353 ieeeok = ilaenv( 10, 'DSTEVR', 'N', 1, 2, 3, 4 )
354*
355 wantz = lsame( jobz, 'V' )
356 alleig = lsame( range, 'A' )
357 valeig = lsame( range, 'V' )
358 indeig = lsame( range, 'I' )
359*
360 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
361 lwmin = max( 1, 20*n )
362 liwmin = max( 1, 10*n )
363*
364*
365 info = 0
366 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
367 info = -1
368 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
369 info = -2
370 ELSE IF( n.LT.0 ) THEN
371 info = -3
372 ELSE
373 IF( valeig ) THEN
374 IF( n.GT.0 .AND. vu.LE.vl )
375 \$ info = -7
376 ELSE IF( indeig ) THEN
377 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
378 info = -8
379 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
380 info = -9
381 END IF
382 END IF
383 END IF
384 IF( info.EQ.0 ) THEN
385 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
386 info = -14
387 END IF
388 END IF
389*
390 IF( info.EQ.0 ) THEN
391 work( 1 ) = lwmin
392 iwork( 1 ) = liwmin
393*
394 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
395 info = -17
396 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
397 info = -19
398 END IF
399 END IF
400*
401 IF( info.NE.0 ) THEN
402 CALL xerbla( 'DSTEVR', -info )
403 RETURN
404 ELSE IF( lquery ) THEN
405 RETURN
406 END IF
407*
408* Quick return if possible
409*
410 m = 0
411 IF( n.EQ.0 )
412 \$ RETURN
413*
414 IF( n.EQ.1 ) THEN
415 IF( alleig .OR. indeig ) THEN
416 m = 1
417 w( 1 ) = d( 1 )
418 ELSE
419 IF( vl.LT.d( 1 ) .AND. vu.GE.d( 1 ) ) THEN
420 m = 1
421 w( 1 ) = d( 1 )
422 END IF
423 END IF
424 IF( wantz )
425 \$ z( 1, 1 ) = one
426 RETURN
427 END IF
428*
429* Get machine constants.
430*
431 safmin = dlamch( 'Safe minimum' )
432 eps = dlamch( 'Precision' )
433 smlnum = safmin / eps
434 bignum = one / smlnum
435 rmin = sqrt( smlnum )
436 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
437*
438*
439* Scale matrix to allowable range, if necessary.
440*
441 iscale = 0
442 IF( valeig ) THEN
443 vll = vl
444 vuu = vu
445 END IF
446*
447 tnrm = dlanst( 'M', n, d, e )
448 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
449 iscale = 1
450 sigma = rmin / tnrm
451 ELSE IF( tnrm.GT.rmax ) THEN
452 iscale = 1
453 sigma = rmax / tnrm
454 END IF
455 IF( iscale.EQ.1 ) THEN
456 CALL dscal( n, sigma, d, 1 )
457 CALL dscal( n-1, sigma, e( 1 ), 1 )
458 IF( valeig ) THEN
459 vll = vl*sigma
460 vuu = vu*sigma
461 END IF
462 END IF
463
464* Initialize indices into workspaces. Note: These indices are used only
465* if DSTERF or DSTEMR fail.
466
467* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
468* stores the block indices of each of the M<=N eigenvalues.
469 indibl = 1
470* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
471* stores the starting and finishing indices of each block.
472 indisp = indibl + n
473* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
474* that corresponding to eigenvectors that fail to converge in
475* DSTEIN. This information is discarded; if any fail, the driver
476* returns INFO > 0.
477 indifl = indisp + n
478* INDIWO is the offset of the remaining integer workspace.
479 indiwo = indisp + n
480*
481* If all eigenvalues are desired, then
482* call DSTERF or DSTEMR. If this fails for some eigenvalue, then
483* try DSTEBZ.
484*
485*
486 test = .false.
487 IF( indeig ) THEN
488 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
489 test = .true.
490 END IF
491 END IF
492 IF( ( alleig .OR. test ) .AND. ieeeok.EQ.1 ) THEN
493 CALL dcopy( n-1, e( 1 ), 1, work( 1 ), 1 )
494 IF( .NOT.wantz ) THEN
495 CALL dcopy( n, d, 1, w, 1 )
496 CALL dsterf( n, w, work, info )
497 ELSE
498 CALL dcopy( n, d, 1, work( n+1 ), 1 )
499 IF (abstol .LE. two*n*eps) THEN
500 tryrac = .true.
501 ELSE
502 tryrac = .false.
503 END IF
504 CALL dstemr( jobz, 'A', n, work( n+1 ), work, vl, vu, il,
505 \$ iu, m, w, z, ldz, n, isuppz, tryrac,
506 \$ work( 2*n+1 ), lwork-2*n, iwork, liwork, info )
507*
508 END IF
509 IF( info.EQ.0 ) THEN
510 m = n
511 GO TO 10
512 END IF
513 info = 0
514 END IF
515*
516* Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN.
517*
518 IF( wantz ) THEN
519 order = 'B'
520 ELSE
521 order = 'E'
522 END IF
523
524 CALL dstebz( range, order, n, vll, vuu, il, iu, abstol, d, e, m,
525 \$ nsplit, w, iwork( indibl ), iwork( indisp ), work,
526 \$ iwork( indiwo ), info )
527*
528 IF( wantz ) THEN
529 CALL dstein( n, d, e, m, w, iwork( indibl ), iwork( indisp ),
530 \$ z, ldz, work, iwork( indiwo ), iwork( indifl ),
531 \$ info )
532 END IF
533*
534* If matrix was scaled, then rescale eigenvalues appropriately.
535*
536 10 CONTINUE
537 IF( iscale.EQ.1 ) THEN
538 IF( info.EQ.0 ) THEN
539 imax = m
540 ELSE
541 imax = info - 1
542 END IF
543 CALL dscal( imax, one / sigma, w, 1 )
544 END IF
545*
546* If eigenvalues are not in order, then sort them, along with
547* eigenvectors.
548*
549 IF( wantz ) THEN
550 DO 30 j = 1, m - 1
551 i = 0
552 tmp1 = w( j )
553 DO 20 jj = j + 1, m
554 IF( w( jj ).LT.tmp1 ) THEN
555 i = jj
556 tmp1 = w( jj )
557 END IF
558 20 CONTINUE
559*
560 IF( i.NE.0 ) THEN
561 itmp1 = iwork( i )
562 w( i ) = w( j )
563 iwork( i ) = iwork( j )
564 w( j ) = tmp1
565 iwork( j ) = itmp1
566 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
567 END IF
568 30 CONTINUE
569 END IF
570*
571* Causes problems with tests 19 & 20:
572* IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002
573*
574*
575 work( 1 ) = lwmin
576 iwork( 1 ) = liwmin
577 RETURN
578*
579* End of DSTEVR
580*
581 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:273
subroutine dstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
DSTEIN
Definition dstein.f:174
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 dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine dstevr(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
DSTEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition dstevr.f:304
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82