LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
ssyevr.f
Go to the documentation of this file.
1*> \brief <b> SSYEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SSYEVR + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyevr.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyevr.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyevr.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
20* ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
21* IWORK, LIWORK, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER JOBZ, RANGE, UPLO
25* INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
26* REAL ABSTOL, VL, VU
27* ..
28* .. Array Arguments ..
29* INTEGER ISUPPZ( * ), IWORK( * )
30* REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SSYEVR computes selected eigenvalues and, optionally, eigenvectors
40*> of a real symmetric matrix A. Eigenvalues and eigenvectors can be
41*> selected by specifying either a range of values or a range of indices
42*> for the desired eigenvalues. Invocations with different choices for
43*> these parameters may result in the computation of slightly different
44*> eigenvalues and/or eigenvectors for the same matrix. The reason for
45*> this behavior is that there exists a variety of algorithms (each
46*> performing best for a particular set of options) with SSYEVR
47*> attempting to select the best based on the various parameters. In all
48*> cases, the computed values are accurate within the limits of finite
49*> precision arithmetic.
50*>
51*> SSYEVR first reduces the matrix A to tridiagonal form T with a call
52*> to SSYTRD. Then, whenever possible, SSYEVR calls SSTEMR to compute
53*> the eigenspectrum using Relatively Robust Representations. SSTEMR
54*> computes eigenvalues by the dqds algorithm, while orthogonal
55*> eigenvectors are computed from various "good" L D L^T representations
56*> (also known as Relatively Robust Representations). Gram-Schmidt
57*> orthogonalization is avoided as far as possible. More specifically,
58*> the various steps of the algorithm are as follows.
59*>
60*> For each unreduced block (submatrix) of T,
61*> (a) Compute T - sigma I = L D L^T, so that L and D
62*> define all the wanted eigenvalues to high relative accuracy.
63*> This means that small relative changes in the entries of D and L
64*> cause only small relative changes in the eigenvalues and
65*> eigenvectors. The standard (unfactored) representation of the
66*> tridiagonal matrix T does not have this property in general.
67*> (b) Compute the eigenvalues to suitable accuracy.
68*> If the eigenvectors are desired, the algorithm attains full
69*> accuracy of the computed eigenvalues only right before
70*> the corresponding vectors have to be computed, see steps c) and d).
71*> (c) For each cluster of close eigenvalues, select a new
72*> shift close to the cluster, find a new factorization, and refine
73*> the shifted eigenvalues to suitable accuracy.
74*> (d) For each eigenvalue with a large enough relative separation compute
75*> the corresponding eigenvector by forming a rank revealing twisted
76*> factorization. Go back to (c) for any clusters that remain.
77*>
78*> The desired accuracy of the output can be specified by the input
79*> parameter ABSTOL.
80*>
81*> For more details, see SSTEMR's documentation and:
82*> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
83*> to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
84*> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
85*> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
86*> Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
87*> 2004. Also LAPACK Working Note 154.
88*> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
89*> tridiagonal eigenvalue/eigenvector problem",
90*> Computer Science Division Technical Report No. UCB/CSD-97-971,
91*> UC Berkeley, May 1997.
92*>
93*>
94*> Note 1 : SSYEVR calls SSTEMR when the full spectrum is requested
95*> on machines which conform to the ieee-754 floating point standard.
96*> SSYEVR calls SSTEBZ and SSTEIN on non-ieee machines and
97*> when partial spectrum requests are made.
98*>
99*> Normal execution of SSTEMR may create NaNs and infinities and
100*> hence may abort due to a floating point exception in environments
101*> which do not handle NaNs and infinities in the ieee standard default
102*> manner.
103*> \endverbatim
104*
105* Arguments:
106* ==========
107*
108*> \param[in] JOBZ
109*> \verbatim
110*> JOBZ is CHARACTER*1
111*> = 'N': Compute eigenvalues only;
112*> = 'V': Compute eigenvalues and eigenvectors.
113*>
114*> This parameter influences the choice of the algorithm and
115*> may alter the computed values.
116*> \endverbatim
117*>
118*> \param[in] RANGE
119*> \verbatim
120*> RANGE is CHARACTER*1
121*> = 'A': all eigenvalues will be found.
122*> = 'V': all eigenvalues in the half-open interval (VL,VU]
123*> will be found.
124*> = 'I': the IL-th through IU-th eigenvalues will be found.
125*> For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
126*> SSTEIN are called
127*>
128*> This parameter influences the choice of the algorithm and
129*> may alter the computed values.
130*> \endverbatim
131*>
132*> \param[in] UPLO
133*> \verbatim
134*> UPLO is CHARACTER*1
135*> = 'U': Upper triangle of A is stored;
136*> = 'L': Lower triangle of A is stored.
137*> \endverbatim
138*>
139*> \param[in] N
140*> \verbatim
141*> N is INTEGER
142*> The order of the matrix A. N >= 0.
143*> \endverbatim
144*>
145*> \param[in,out] A
146*> \verbatim
147*> A is REAL array, dimension (LDA, N)
148*> On entry, the symmetric matrix A. If UPLO = 'U', the
149*> leading N-by-N upper triangular part of A contains the
150*> upper triangular part of the matrix A. If UPLO = 'L',
151*> the leading N-by-N lower triangular part of A contains
152*> the lower triangular part of the matrix A.
153*> On exit, the lower triangle (if UPLO='L') or the upper
154*> triangle (if UPLO='U') of A, including the diagonal, is
155*> destroyed.
156*> \endverbatim
157*>
158*> \param[in] LDA
159*> \verbatim
160*> LDA is INTEGER
161*> The leading dimension of the array A. LDA >= max(1,N).
162*> \endverbatim
163*>
164*> \param[in] VL
165*> \verbatim
166*> VL is REAL
167*> If RANGE='V', the lower bound of the interval to
168*> be searched for eigenvalues. VL < VU.
169*> Not referenced if RANGE = 'A' or 'I'.
170*> \endverbatim
171*>
172*> \param[in] VU
173*> \verbatim
174*> VU is REAL
175*> If RANGE='V', the upper bound of the interval to
176*> be searched for eigenvalues. VL < VU.
177*> Not referenced if RANGE = 'A' or 'I'.
178*> \endverbatim
179*>
180*> \param[in] IL
181*> \verbatim
182*> IL is INTEGER
183*> If RANGE='I', the index of the
184*> smallest eigenvalue to be returned.
185*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
186*> Not referenced if RANGE = 'A' or 'V'.
187*> \endverbatim
188*>
189*> \param[in] IU
190*> \verbatim
191*> IU is INTEGER
192*> If RANGE='I', the index of the
193*> largest eigenvalue to be returned.
194*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
195*> Not referenced if RANGE = 'A' or 'V'.
196*> \endverbatim
197*>
198*> \param[in] ABSTOL
199*> \verbatim
200*> ABSTOL is REAL
201*> The absolute error tolerance for the eigenvalues.
202*> An approximate eigenvalue is accepted as converged
203*> when it is determined to lie in an interval [a,b]
204*> of width less than or equal to
205*>
206*> ABSTOL + EPS * max( |a|,|b| ) ,
207*>
208*> where EPS is the machine precision. If ABSTOL is less than
209*> or equal to zero, then EPS*|T| will be used in its place,
210*> where |T| is the 1-norm of the tridiagonal matrix obtained
211*> by reducing A to tridiagonal form.
212*>
213*> See "Computing Small Singular Values of Bidiagonal Matrices
214*> with Guaranteed High Relative Accuracy," by Demmel and
215*> Kahan, LAPACK Working Note #3.
216*>
217*> If high relative accuracy is important, set ABSTOL to
218*> SLAMCH( 'Safe minimum' ). Doing so will guarantee that
219*> eigenvalues are computed to high relative accuracy when
220*> possible in future releases. The current code does not
221*> make any guarantees about high relative accuracy, but
222*> future releases will. See J. Barlow and J. Demmel,
223*> "Computing Accurate Eigensystems of Scaled Diagonally
224*> Dominant Matrices", LAPACK Working Note #7, for a discussion
225*> of which matrices define their eigenvalues to high relative
226*> accuracy.
227*> \endverbatim
228*>
229*> \param[out] M
230*> \verbatim
231*> M is INTEGER
232*> The total number of eigenvalues found. 0 <= M <= N.
233*> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
234*> \endverbatim
235*>
236*> \param[out] W
237*> \verbatim
238*> W is REAL array, dimension (N)
239*> The first M elements contain the selected eigenvalues in
240*> ascending order.
241*> \endverbatim
242*>
243*> \param[out] Z
244*> \verbatim
245*> Z is REAL array, dimension (LDZ, max(1,M))
246*> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
247*> contain the orthonormal eigenvectors of the matrix A
248*> corresponding to the selected eigenvalues, with the i-th
249*> column of Z holding the eigenvector associated with W(i).
250*> If JOBZ = 'N', then Z is not referenced.
251*> Note: the user must ensure that at least max(1,M) columns are
252*> supplied in the array Z; if RANGE = 'V', the exact value of M
253*> is not known in advance and an upper bound must be used.
254*> Supplying N columns is always safe.
255*> \endverbatim
256*>
257*> \param[in] LDZ
258*> \verbatim
259*> LDZ is INTEGER
260*> The leading dimension of the array Z. LDZ >= 1, and if
261*> JOBZ = 'V', LDZ >= max(1,N).
262*> \endverbatim
263*>
264*> \param[out] ISUPPZ
265*> \verbatim
266*> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
267*> The support of the eigenvectors in Z, i.e., the indices
268*> indicating the nonzero elements in Z. The i-th eigenvector
269*> is nonzero only in elements ISUPPZ( 2*i-1 ) through
270*> ISUPPZ( 2*i ). This is an output of SSTEMR (tridiagonal
271*> matrix). The support of the eigenvectors of A is typically
272*> 1:N because of the orthogonal transformations applied by SORMTR.
273*> Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
274*> \endverbatim
275*>
276*> \param[out] WORK
277*> \verbatim
278*> WORK is REAL array, dimension (MAX(1,LWORK))
279*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
280*> \endverbatim
281*>
282*> \param[in] LWORK
283*> \verbatim
284*> LWORK is INTEGER
285*> The dimension of the array WORK.
286*> If N <= 1, LWORK >= 1, else LWORK >= 26*N.
287*> For optimal efficiency, LWORK >= (NB+6)*N,
288*> where NB is the max of the blocksize for SSYTRD and SORMTR
289*> returned by ILAENV.
290*>
291*> If LWORK = -1, then a workspace query is assumed; the routine
292*> only calculates the optimal sizes of the WORK and IWORK
293*> arrays, returns these values as the first entries of the WORK
294*> and IWORK arrays, and no error message related to LWORK or
295*> LIWORK is issued by XERBLA.
296*> \endverbatim
297*>
298*> \param[out] IWORK
299*> \verbatim
300*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
301*> On exit, if INFO = 0, IWORK(1) returns the optimal LWORK.
302*> \endverbatim
303*>
304*> \param[in] LIWORK
305*> \verbatim
306*> LIWORK is INTEGER
307*> The dimension of the array IWORK.
308*> If N <= 1, LIWORK >= 1, else LIWORK >= 10*N.
309*>
310*> If LIWORK = -1, then a workspace query is assumed; the
311*> routine only calculates the optimal sizes of the WORK and
312*> IWORK arrays, returns these values as the first entries of
313*> the WORK and IWORK arrays, and no error message related to
314*> LWORK or LIWORK is issued by XERBLA.
315*> \endverbatim
316*>
317*> \param[out] INFO
318*> \verbatim
319*> INFO is INTEGER
320*> = 0: successful exit
321*> < 0: if INFO = -i, the i-th argument had an illegal value
322*> > 0: Internal error
323*> \endverbatim
324*
325* Authors:
326* ========
327*
328*> \author Univ. of Tennessee
329*> \author Univ. of California Berkeley
330*> \author Univ. of Colorado Denver
331*> \author NAG Ltd.
332*
333*> \ingroup heevr
334*
335*> \par Contributors:
336* ==================
337*>
338*> Inderjit Dhillon, IBM Almaden, USA \n
339*> Osni Marques, LBNL/NERSC, USA \n
340*> Ken Stanley, Computer Science Division, University of
341*> California at Berkeley, USA \n
342*> Jason Riedy, Computer Science Division, University of
343*> California at Berkeley, USA \n
344*>
345* =====================================================================
346 SUBROUTINE ssyevr( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL,
347 $ IU,
348 $ ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
349 $ IWORK, LIWORK, INFO )
350*
351* -- LAPACK driver routine --
352* -- LAPACK is a software package provided by Univ. of Tennessee, --
353* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
354*
355* .. Scalar Arguments ..
356 CHARACTER JOBZ, RANGE, UPLO
357 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
358 REAL ABSTOL, VL, VU
359* ..
360* .. Array Arguments ..
361 INTEGER ISUPPZ( * ), IWORK( * )
362 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
363* ..
364*
365* =====================================================================
366*
367* .. Parameters ..
368 REAL ZERO, ONE, TWO
369 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
370* ..
371* .. Local Scalars ..
372 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
373 $ WANTZ, TRYRAC
374 CHARACTER ORDER
375 INTEGER I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
376 $ INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
377 $ indwk, indwkn, iscale, j, jj, liwmin,
378 $ llwork, llwrkn, lwkopt, lwmin, nb, nsplit
379 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
380 $ SIGMA, SMLNUM, TMP1, VLL, VUU
381* ..
382* .. External Functions ..
383 LOGICAL LSAME
384 INTEGER ILAENV
385 REAL SLAMCH, SLANSY, SROUNDUP_LWORK
386 EXTERNAL lsame, ilaenv, slamch,
387 $ slansy, sroundup_lwork
388* ..
389* .. External Subroutines ..
390 EXTERNAL scopy, sormtr, sscal, sstebz, sstemr,
391 $ sstein,
393* ..
394* .. Intrinsic Functions ..
395 INTRINSIC max, min, sqrt
396* ..
397* .. Executable Statements ..
398*
399* Test the input parameters.
400*
401 ieeeok = ilaenv( 10, 'SSYEVR', 'N', 1, 2, 3, 4 )
402*
403 lower = lsame( uplo, 'L' )
404 wantz = lsame( jobz, 'V' )
405 alleig = lsame( range, 'A' )
406 valeig = lsame( range, 'V' )
407 indeig = lsame( range, 'I' )
408*
409 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
410*
411 IF( n.LE.1 ) THEN
412 lwmin = 1
413 liwmin = 1
414 ELSE
415 lwmin = 26*n
416 liwmin = 10*n
417 END IF
418*
419 info = 0
420 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
421 info = -1
422 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
423 info = -2
424 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
425 info = -3
426 ELSE IF( n.LT.0 ) THEN
427 info = -4
428 ELSE IF( lda.LT.max( 1, n ) ) THEN
429 info = -6
430 ELSE
431 IF( valeig ) THEN
432 IF( n.GT.0 .AND. vu.LE.vl )
433 $ info = -8
434 ELSE IF( indeig ) THEN
435 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
436 info = -9
437 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
438 info = -10
439 END IF
440 END IF
441 END IF
442 IF( info.EQ.0 ) THEN
443 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
444 info = -15
445 END IF
446 END IF
447*
448 IF( info.EQ.0 ) THEN
449 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
450 nb = max( nb, ilaenv( 1, 'SORMTR', uplo, n, -1, -1, -1 ) )
451 lwkopt = max( ( nb+1 )*n, lwmin )
452 work( 1 ) = sroundup_lwork( lwkopt )
453 iwork( 1 ) = liwmin
454*
455 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
456 info = -18
457 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
458 info = -20
459 END IF
460 END IF
461*
462 IF( info.NE.0 ) THEN
463 CALL xerbla( 'SSYEVR', -info )
464 RETURN
465 ELSE IF( lquery ) THEN
466 RETURN
467 END IF
468*
469* Quick return if possible
470*
471 m = 0
472 IF( n.EQ.0 ) THEN
473 work( 1 ) = 1
474 RETURN
475 END IF
476*
477 IF( n.EQ.1 ) THEN
478 work( 1 ) = 26
479 IF( alleig .OR. indeig ) THEN
480 m = 1
481 w( 1 ) = a( 1, 1 )
482 ELSE
483 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
484 m = 1
485 w( 1 ) = a( 1, 1 )
486 END IF
487 END IF
488 IF( wantz ) THEN
489 z( 1, 1 ) = one
490 isuppz( 1 ) = 1
491 isuppz( 2 ) = 1
492 END IF
493 RETURN
494 END IF
495*
496* Get machine constants.
497*
498 safmin = slamch( 'Safe minimum' )
499 eps = slamch( 'Precision' )
500 smlnum = safmin / eps
501 bignum = one / smlnum
502 rmin = sqrt( smlnum )
503 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
504*
505* Scale matrix to allowable range, if necessary.
506*
507 iscale = 0
508 abstll = abstol
509 IF (valeig) THEN
510 vll = vl
511 vuu = vu
512 END IF
513 anrm = slansy( 'M', uplo, n, a, lda, work )
514 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
515 iscale = 1
516 sigma = rmin / anrm
517 ELSE IF( anrm.GT.rmax ) THEN
518 iscale = 1
519 sigma = rmax / anrm
520 END IF
521 IF( iscale.EQ.1 ) THEN
522 IF( lower ) THEN
523 DO 10 j = 1, n
524 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
525 10 CONTINUE
526 ELSE
527 DO 20 j = 1, n
528 CALL sscal( j, sigma, a( 1, j ), 1 )
529 20 CONTINUE
530 END IF
531 IF( abstol.GT.0 )
532 $ abstll = abstol*sigma
533 IF( valeig ) THEN
534 vll = vl*sigma
535 vuu = vu*sigma
536 END IF
537 END IF
538
539* Initialize indices into workspaces. Note: The IWORK indices are
540* used only if SSTERF or SSTEMR fail.
541
542* WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
543* elementary reflectors used in SSYTRD.
544 indtau = 1
545* WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
546 indd = indtau + n
547* WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
548* tridiagonal matrix from SSYTRD.
549 inde = indd + n
550* WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
551* -written by SSTEMR (the SSTERF path copies the diagonal to W).
552 inddd = inde + n
553* WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
554* -written while computing the eigenvalues in SSTERF and SSTEMR.
555 indee = inddd + n
556* INDWK is the starting offset of the left-over workspace, and
557* LLWORK is the remaining workspace size.
558 indwk = indee + n
559 llwork = lwork - indwk + 1
560
561* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
562* stores the block indices of each of the M<=N eigenvalues.
563 indibl = 1
564* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
565* stores the starting and finishing indices of each block.
566 indisp = indibl + n
567* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
568* that corresponding to eigenvectors that fail to converge in
569* SSTEIN. This information is discarded; if any fail, the driver
570* returns INFO > 0.
571 indifl = indisp + n
572* INDIWO is the offset of the remaining integer workspace.
573 indiwo = indifl + n
574
575*
576* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
577*
578 CALL ssytrd( uplo, n, a, lda, work( indd ), work( inde ),
579 $ work( indtau ), work( indwk ), llwork, iinfo )
580*
581* If all eigenvalues are desired
582* then call SSTERF or SSTEMR and SORMTR.
583*
584 test = .false.
585 IF( indeig ) THEN
586 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
587 test = .true.
588 END IF
589 END IF
590 IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
591 IF( .NOT.wantz ) THEN
592 CALL scopy( n, work( indd ), 1, w, 1 )
593 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
594 CALL ssterf( n, w, work( indee ), info )
595 ELSE
596 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
597 CALL scopy( n, work( indd ), 1, work( inddd ), 1 )
598*
599 IF (abstol .LE. two*real( n )*eps) THEN
600 tryrac = .true.
601 ELSE
602 tryrac = .false.
603 END IF
604 CALL sstemr( jobz, 'A', n, work( inddd ), work( indee ),
605 $ vl, vu, il, iu, m, w, z, ldz, n, isuppz,
606 $ tryrac, work( indwk ), lwork, iwork, liwork,
607 $ info )
608*
609*
610*
611* Apply orthogonal matrix used in reduction to tridiagonal
612* form to eigenvectors returned by SSTEMR.
613*
614 IF( wantz .AND. info.EQ.0 ) THEN
615 indwkn = inde
616 llwrkn = lwork - indwkn + 1
617 CALL sormtr( 'L', uplo, 'N', n, m, a, lda,
618 $ work( indtau ), z, ldz, work( indwkn ),
619 $ llwrkn, iinfo )
620 END IF
621 END IF
622*
623*
624 IF( info.EQ.0 ) THEN
625* Everything worked. Skip SSTEBZ/SSTEIN. IWORK(:) are
626* undefined.
627 m = n
628 GO TO 30
629 END IF
630 info = 0
631 END IF
632*
633* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
634* Also call SSTEBZ and SSTEIN if SSTEMR fails.
635*
636 IF( wantz ) THEN
637 order = 'B'
638 ELSE
639 order = 'E'
640 END IF
641
642 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
643 $ work( indd ), work( inde ), m, nsplit, w,
644 $ iwork( indibl ), iwork( indisp ), work( indwk ),
645 $ iwork( indiwo ), info )
646*
647 IF( wantz ) THEN
648 CALL sstein( n, work( indd ), work( inde ), m, w,
649 $ iwork( indibl ), iwork( indisp ), z, ldz,
650 $ work( indwk ), iwork( indiwo ), iwork( indifl ),
651 $ info )
652*
653* Apply orthogonal matrix used in reduction to tridiagonal
654* form to eigenvectors returned by SSTEIN.
655*
656 indwkn = inde
657 llwrkn = lwork - indwkn + 1
658 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ),
659 $ z,
660 $ ldz, work( indwkn ), llwrkn, iinfo )
661 END IF
662*
663* If matrix was scaled, then rescale eigenvalues appropriately.
664*
665* Jump here if SSTEMR/SSTEIN succeeded.
666 30 CONTINUE
667 IF( iscale.EQ.1 ) THEN
668 IF( info.EQ.0 ) THEN
669 imax = m
670 ELSE
671 imax = info - 1
672 END IF
673 CALL sscal( imax, one / sigma, w, 1 )
674 END IF
675*
676* If eigenvalues are not in order, then sort them, along with
677* eigenvectors. Note: We do not sort the IFAIL portion of IWORK.
678* It may not be initialized (if SSTEMR/SSTEIN succeeded), and we do
679* not return this detailed information to the user.
680*
681 IF( wantz ) THEN
682 DO 50 j = 1, m - 1
683 i = 0
684 tmp1 = w( j )
685 DO 40 jj = j + 1, m
686 IF( w( jj ).LT.tmp1 ) THEN
687 i = jj
688 tmp1 = w( jj )
689 END IF
690 40 CONTINUE
691*
692 IF( i.NE.0 ) THEN
693 w( i ) = w( j )
694 w( j ) = tmp1
695 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
696 END IF
697 50 CONTINUE
698 END IF
699*
700* Set WORK(1) to optimal workspace size.
701*
702 work( 1 ) = sroundup_lwork( lwkopt )
703 iwork( 1 ) = liwmin
704*
705 RETURN
706*
707* End of SSYEVR
708*
709 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ssyevr(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
SSYEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition ssyevr.f:350
subroutine ssytrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
SSYTRD
Definition ssytrd.f:191
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:272
subroutine sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:172
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:320
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:84
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
subroutine sormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
SORMTR
Definition sormtr.f:171