LAPACK  3.7.1
LAPACK: Linear Algebra PACKage
ssyevx.f
Go to the documentation of this file.
1 *> \brief <b> SSYEVX 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 *> \htmlonly
9 *> Download SSYEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
22 * ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK,
23 * IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
28 * REAL ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> SSYEVX computes selected eigenvalues and, optionally, eigenvectors
42 *> of a real symmetric matrix A. Eigenvalues and eigenvectors can be
43 *> selected by specifying either a range of values or a range of indices
44 *> for the desired eigenvalues.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] JOBZ
51 *> \verbatim
52 *> JOBZ is CHARACTER*1
53 *> = 'N': Compute eigenvalues only;
54 *> = 'V': Compute eigenvalues and eigenvectors.
55 *> \endverbatim
56 *>
57 *> \param[in] RANGE
58 *> \verbatim
59 *> RANGE is CHARACTER*1
60 *> = 'A': all eigenvalues will be found.
61 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
62 *> will be found.
63 *> = 'I': the IL-th through IU-th eigenvalues will be found.
64 *> \endverbatim
65 *>
66 *> \param[in] UPLO
67 *> \verbatim
68 *> UPLO is CHARACTER*1
69 *> = 'U': Upper triangle of A is stored;
70 *> = 'L': Lower triangle of A is stored.
71 *> \endverbatim
72 *>
73 *> \param[in] N
74 *> \verbatim
75 *> N is INTEGER
76 *> The order of the matrix A. N >= 0.
77 *> \endverbatim
78 *>
79 *> \param[in,out] A
80 *> \verbatim
81 *> A is REAL array, dimension (LDA, N)
82 *> On entry, the symmetric matrix A. If UPLO = 'U', the
83 *> leading N-by-N upper triangular part of A contains the
84 *> upper triangular part of the matrix A. If UPLO = 'L',
85 *> the leading N-by-N lower triangular part of A contains
86 *> the lower triangular part of the matrix A.
87 *> On exit, the lower triangle (if UPLO='L') or the upper
88 *> triangle (if UPLO='U') of A, including the diagonal, is
89 *> destroyed.
90 *> \endverbatim
91 *>
92 *> \param[in] LDA
93 *> \verbatim
94 *> LDA is INTEGER
95 *> The leading dimension of the array A. LDA >= max(1,N).
96 *> \endverbatim
97 *>
98 *> \param[in] VL
99 *> \verbatim
100 *> VL is REAL
101 *> If RANGE='V', the lower bound of the interval to
102 *> be searched for eigenvalues. VL < VU.
103 *> Not referenced if RANGE = 'A' or 'I'.
104 *> \endverbatim
105 *>
106 *> \param[in] VU
107 *> \verbatim
108 *> VU is REAL
109 *> If RANGE='V', the upper bound of the interval to
110 *> be searched for eigenvalues. VL < VU.
111 *> Not referenced if RANGE = 'A' or 'I'.
112 *> \endverbatim
113 *>
114 *> \param[in] IL
115 *> \verbatim
116 *> IL is INTEGER
117 *> If RANGE='I', the index of the
118 *> smallest eigenvalue to be returned.
119 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
120 *> Not referenced if RANGE = 'A' or 'V'.
121 *> \endverbatim
122 *>
123 *> \param[in] IU
124 *> \verbatim
125 *> IU is INTEGER
126 *> If RANGE='I', the index of the
127 *> largest eigenvalue to be returned.
128 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
129 *> Not referenced if RANGE = 'A' or 'V'.
130 *> \endverbatim
131 *>
132 *> \param[in] ABSTOL
133 *> \verbatim
134 *> ABSTOL is REAL
135 *> The absolute error tolerance for the eigenvalues.
136 *> An approximate eigenvalue is accepted as converged
137 *> when it is determined to lie in an interval [a,b]
138 *> of width less than or equal to
139 *>
140 *> ABSTOL + EPS * max( |a|,|b| ) ,
141 *>
142 *> where EPS is the machine precision. If ABSTOL is less than
143 *> or equal to zero, then EPS*|T| will be used in its place,
144 *> where |T| is the 1-norm of the tridiagonal matrix obtained
145 *> by reducing A to tridiagonal form.
146 *>
147 *> Eigenvalues will be computed most accurately when ABSTOL is
148 *> set to twice the underflow threshold 2*SLAMCH('S'), not zero.
149 *> If this routine returns with INFO>0, indicating that some
150 *> eigenvectors did not converge, try setting ABSTOL to
151 *> 2*SLAMCH('S').
152 *>
153 *> See "Computing Small Singular Values of Bidiagonal Matrices
154 *> with Guaranteed High Relative Accuracy," by Demmel and
155 *> Kahan, LAPACK Working Note #3.
156 *> \endverbatim
157 *>
158 *> \param[out] M
159 *> \verbatim
160 *> M is INTEGER
161 *> The total number of eigenvalues found. 0 <= M <= N.
162 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
163 *> \endverbatim
164 *>
165 *> \param[out] W
166 *> \verbatim
167 *> W is REAL array, dimension (N)
168 *> On normal exit, the first M elements contain the selected
169 *> eigenvalues in ascending order.
170 *> \endverbatim
171 *>
172 *> \param[out] Z
173 *> \verbatim
174 *> Z is REAL array, dimension (LDZ, max(1,M))
175 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
176 *> contain the orthonormal eigenvectors of the matrix A
177 *> corresponding to the selected eigenvalues, with the i-th
178 *> column of Z holding the eigenvector associated with W(i).
179 *> If an eigenvector fails to converge, then that column of Z
180 *> contains the latest approximation to the eigenvector, and the
181 *> index of the eigenvector is returned in IFAIL.
182 *> If JOBZ = 'N', then Z is not referenced.
183 *> Note: the user must ensure that at least max(1,M) columns are
184 *> supplied in the array Z; if RANGE = 'V', the exact value of M
185 *> is not known in advance and an upper bound must be used.
186 *> \endverbatim
187 *>
188 *> \param[in] LDZ
189 *> \verbatim
190 *> LDZ is INTEGER
191 *> The leading dimension of the array Z. LDZ >= 1, and if
192 *> JOBZ = 'V', LDZ >= max(1,N).
193 *> \endverbatim
194 *>
195 *> \param[out] WORK
196 *> \verbatim
197 *> WORK is REAL array, dimension (MAX(1,LWORK))
198 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
199 *> \endverbatim
200 *>
201 *> \param[in] LWORK
202 *> \verbatim
203 *> LWORK is INTEGER
204 *> The length of the array WORK. LWORK >= 1, when N <= 1;
205 *> otherwise 8*N.
206 *> For optimal efficiency, LWORK >= (NB+3)*N,
207 *> where NB is the max of the blocksize for SSYTRD and SORMTR
208 *> returned by ILAENV.
209 *>
210 *> If LWORK = -1, then a workspace query is assumed; the routine
211 *> only calculates the optimal size of the WORK array, returns
212 *> this value as the first entry of the WORK array, and no error
213 *> message related to LWORK is issued by XERBLA.
214 *> \endverbatim
215 *>
216 *> \param[out] IWORK
217 *> \verbatim
218 *> IWORK is INTEGER array, dimension (5*N)
219 *> \endverbatim
220 *>
221 *> \param[out] IFAIL
222 *> \verbatim
223 *> IFAIL is INTEGER array, dimension (N)
224 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
225 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
226 *> indices of the eigenvectors that failed to converge.
227 *> If JOBZ = 'N', then IFAIL is not referenced.
228 *> \endverbatim
229 *>
230 *> \param[out] INFO
231 *> \verbatim
232 *> INFO is INTEGER
233 *> = 0: successful exit
234 *> < 0: if INFO = -i, the i-th argument had an illegal value
235 *> > 0: if INFO = i, then i eigenvectors failed to converge.
236 *> Their indices are stored in array IFAIL.
237 *> \endverbatim
238 *
239 * Authors:
240 * ========
241 *
242 *> \author Univ. of Tennessee
243 *> \author Univ. of California Berkeley
244 *> \author Univ. of Colorado Denver
245 *> \author NAG Ltd.
246 *
247 *> \date June 2016
248 *
249 *> \ingroup realSYeigen
250 *
251 * =====================================================================
252  SUBROUTINE ssyevx( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
253  $ ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK,
254  $ IFAIL, INFO )
255 *
256 * -- LAPACK driver routine (version 3.7.0) --
257 * -- LAPACK is a software package provided by Univ. of Tennessee, --
258 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
259 * June 2016
260 *
261 * .. Scalar Arguments ..
262  CHARACTER JOBZ, RANGE, UPLO
263  INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
264  REAL ABSTOL, VL, VU
265 * ..
266 * .. Array Arguments ..
267  INTEGER IFAIL( * ), IWORK( * )
268  REAL A( lda, * ), W( * ), WORK( * ), Z( ldz, * )
269 * ..
270 *
271 * =====================================================================
272 *
273 * .. Parameters ..
274  REAL ZERO, ONE
275  parameter( zero = 0.0e+0, one = 1.0e+0 )
276 * ..
277 * .. Local Scalars ..
278  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
279  $ wantz
280  CHARACTER ORDER
281  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
282  $ indisp, indiwo, indtau, indwkn, indwrk, iscale,
283  $ itmp1, j, jj, llwork, llwrkn, lwkmin,
284  $ lwkopt, nb, nsplit
285  REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
286  $ sigma, smlnum, tmp1, vll, vuu
287 * ..
288 * .. External Functions ..
289  LOGICAL LSAME
290  INTEGER ILAENV
291  REAL SLAMCH, SLANSY
292  EXTERNAL lsame, ilaenv, slamch, slansy
293 * ..
294 * .. External Subroutines ..
295  EXTERNAL scopy, slacpy, sorgtr, sormtr, sscal, sstebz,
297 * ..
298 * .. Intrinsic Functions ..
299  INTRINSIC max, min, sqrt
300 * ..
301 * .. Executable Statements ..
302 *
303 * Test the input parameters.
304 *
305  lower = lsame( uplo, 'L' )
306  wantz = lsame( jobz, 'V' )
307  alleig = lsame( range, 'A' )
308  valeig = lsame( range, 'V' )
309  indeig = lsame( range, 'I' )
310  lquery = ( lwork.EQ.-1 )
311 *
312  info = 0
313  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
314  info = -1
315  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
316  info = -2
317  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
318  info = -3
319  ELSE IF( n.LT.0 ) THEN
320  info = -4
321  ELSE IF( lda.LT.max( 1, n ) ) THEN
322  info = -6
323  ELSE
324  IF( valeig ) THEN
325  IF( n.GT.0 .AND. vu.LE.vl )
326  $ info = -8
327  ELSE IF( indeig ) THEN
328  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
329  info = -9
330  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
331  info = -10
332  END IF
333  END IF
334  END IF
335  IF( info.EQ.0 ) THEN
336  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
337  info = -15
338  END IF
339  END IF
340 *
341  IF( info.EQ.0 ) THEN
342  IF( n.LE.1 ) THEN
343  lwkmin = 1
344  work( 1 ) = lwkmin
345  ELSE
346  lwkmin = 8*n
347  nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
348  nb = max( nb, ilaenv( 1, 'SORMTR', uplo, n, -1, -1, -1 ) )
349  lwkopt = max( lwkmin, ( nb + 3 )*n )
350  work( 1 ) = lwkopt
351  END IF
352 *
353  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
354  $ info = -17
355  END IF
356 *
357  IF( info.NE.0 ) THEN
358  CALL xerbla( 'SSYEVX', -info )
359  RETURN
360  ELSE IF( lquery ) THEN
361  RETURN
362  END IF
363 *
364 * Quick return if possible
365 *
366  m = 0
367  IF( n.EQ.0 ) THEN
368  RETURN
369  END IF
370 *
371  IF( n.EQ.1 ) THEN
372  IF( alleig .OR. indeig ) THEN
373  m = 1
374  w( 1 ) = a( 1, 1 )
375  ELSE
376  IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
377  m = 1
378  w( 1 ) = a( 1, 1 )
379  END IF
380  END IF
381  IF( wantz )
382  $ z( 1, 1 ) = one
383  RETURN
384  END IF
385 *
386 * Get machine constants.
387 *
388  safmin = slamch( 'Safe minimum' )
389  eps = slamch( 'Precision' )
390  smlnum = safmin / eps
391  bignum = one / smlnum
392  rmin = sqrt( smlnum )
393  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
394 *
395 * Scale matrix to allowable range, if necessary.
396 *
397  iscale = 0
398  abstll = abstol
399  IF( valeig ) THEN
400  vll = vl
401  vuu = vu
402  END IF
403  anrm = slansy( 'M', uplo, n, a, lda, work )
404  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
405  iscale = 1
406  sigma = rmin / anrm
407  ELSE IF( anrm.GT.rmax ) THEN
408  iscale = 1
409  sigma = rmax / anrm
410  END IF
411  IF( iscale.EQ.1 ) THEN
412  IF( lower ) THEN
413  DO 10 j = 1, n
414  CALL sscal( n-j+1, sigma, a( j, j ), 1 )
415  10 CONTINUE
416  ELSE
417  DO 20 j = 1, n
418  CALL sscal( j, sigma, a( 1, j ), 1 )
419  20 CONTINUE
420  END IF
421  IF( abstol.GT.0 )
422  $ abstll = abstol*sigma
423  IF( valeig ) THEN
424  vll = vl*sigma
425  vuu = vu*sigma
426  END IF
427  END IF
428 *
429 * Call SSYTRD to reduce symmetric matrix to tridiagonal form.
430 *
431  indtau = 1
432  inde = indtau + n
433  indd = inde + n
434  indwrk = indd + n
435  llwork = lwork - indwrk + 1
436  CALL ssytrd( uplo, n, a, lda, work( indd ), work( inde ),
437  $ work( indtau ), work( indwrk ), llwork, iinfo )
438 *
439 * If all eigenvalues are desired and ABSTOL is less than or equal to
440 * zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
441 * some eigenvalue, then try SSTEBZ.
442 *
443  test = .false.
444  IF( indeig ) THEN
445  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
446  test = .true.
447  END IF
448  END IF
449  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
450  CALL scopy( n, work( indd ), 1, w, 1 )
451  indee = indwrk + 2*n
452  IF( .NOT.wantz ) THEN
453  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
454  CALL ssterf( n, w, work( indee ), info )
455  ELSE
456  CALL slacpy( 'A', n, n, a, lda, z, ldz )
457  CALL sorgtr( uplo, n, z, ldz, work( indtau ),
458  $ work( indwrk ), llwork, iinfo )
459  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
460  CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
461  $ work( indwrk ), info )
462  IF( info.EQ.0 ) THEN
463  DO 30 i = 1, n
464  ifail( i ) = 0
465  30 CONTINUE
466  END IF
467  END IF
468  IF( info.EQ.0 ) THEN
469  m = n
470  GO TO 40
471  END IF
472  info = 0
473  END IF
474 *
475 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
476 *
477  IF( wantz ) THEN
478  order = 'B'
479  ELSE
480  order = 'E'
481  END IF
482  indibl = 1
483  indisp = indibl + n
484  indiwo = indisp + n
485  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
486  $ work( indd ), work( inde ), m, nsplit, w,
487  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
488  $ iwork( indiwo ), info )
489 *
490  IF( wantz ) THEN
491  CALL sstein( n, work( indd ), work( inde ), m, w,
492  $ iwork( indibl ), iwork( indisp ), z, ldz,
493  $ work( indwrk ), iwork( indiwo ), ifail, info )
494 *
495 * Apply orthogonal matrix used in reduction to tridiagonal
496 * form to eigenvectors returned by SSTEIN.
497 *
498  indwkn = inde
499  llwrkn = lwork - indwkn + 1
500  CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
501  $ ldz, work( indwkn ), llwrkn, iinfo )
502  END IF
503 *
504 * If matrix was scaled, then rescale eigenvalues appropriately.
505 *
506  40 CONTINUE
507  IF( iscale.EQ.1 ) THEN
508  IF( info.EQ.0 ) THEN
509  imax = m
510  ELSE
511  imax = info - 1
512  END IF
513  CALL sscal( imax, one / sigma, w, 1 )
514  END IF
515 *
516 * If eigenvalues are not in order, then sort them, along with
517 * eigenvectors.
518 *
519  IF( wantz ) THEN
520  DO 60 j = 1, m - 1
521  i = 0
522  tmp1 = w( j )
523  DO 50 jj = j + 1, m
524  IF( w( jj ).LT.tmp1 ) THEN
525  i = jj
526  tmp1 = w( jj )
527  END IF
528  50 CONTINUE
529 *
530  IF( i.NE.0 ) THEN
531  itmp1 = iwork( indibl+i-1 )
532  w( i ) = w( j )
533  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
534  w( j ) = tmp1
535  iwork( indibl+j-1 ) = itmp1
536  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
537  IF( info.NE.0 ) THEN
538  itmp1 = ifail( i )
539  ifail( i ) = ifail( j )
540  ifail( j ) = itmp1
541  END IF
542  END IF
543  60 CONTINUE
544  END IF
545 *
546 * Set WORK(1) to optimal workspace size.
547 *
548  work( 1 ) = lwkopt
549 *
550  RETURN
551 *
552 * End of SSYEVX
553 *
554  END
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:275
subroutine ssyevx(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO)
SSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices ...
Definition: ssyevx.f:255
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine ssytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
SSYTRD
Definition: ssytrd.f:194
subroutine sormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMTR
Definition: sormtr.f:174
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:176
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine sorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
SORGTR
Definition: sorgtr.f:125
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84