LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ssyevx()

subroutine ssyevx ( character  jobz,
character  range,
character  uplo,
integer  n,
real, dimension( lda, * )  a,
integer  lda,
real  vl,
real  vu,
integer  il,
integer  iu,
real  abstol,
integer  m,
real, dimension( * )  w,
real, dimension( ldz, * )  z,
integer  ldz,
real, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer, dimension( * )  ifail,
integer  info 
)

SSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

Download SSYEVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SSYEVX computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
 selected by specifying either a range of values or a range of indices
 for the desired eigenvalues.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found.
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th eigenvalues will be found.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is REAL array, dimension (LDA, N)
          On entry, the symmetric matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]VL
          VL is REAL
          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
          VU is REAL
          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]ABSTOL
          ABSTOL is REAL
          The absolute error tolerance for the eigenvalues.
          An approximate eigenvalue is accepted as converged
          when it is determined to lie in an interval [a,b]
          of width less than or equal to

                  ABSTOL + EPS *   max( |a|,|b| ) ,

          where EPS is the machine precision.  If ABSTOL is less than
          or equal to zero, then  EPS*|T|  will be used in its place,
          where |T| is the 1-norm of the tridiagonal matrix obtained
          by reducing A to tridiagonal form.

          Eigenvalues will be computed most accurately when ABSTOL is
          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*SLAMCH('S').

          See "Computing Small Singular Values of Bidiagonal Matrices
          with Guaranteed High Relative Accuracy," by Demmel and
          Kahan, LAPACK Working Note #3.
[out]M
          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
[out]W
          W is REAL array, dimension (N)
          On normal exit, the first M elements contain the selected
          eigenvalues in ascending order.
[out]Z
          Z is REAL array, dimension (LDZ, max(1,M))
          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
          contain the orthonormal eigenvectors of the matrix A
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          If an eigenvector fails to converge, then that column of Z
          contains the latest approximation to the eigenvector, and the
          index of the eigenvector is returned in IFAIL.
          If JOBZ = 'N', then Z is not referenced.
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and an upper bound must be used.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 1, when N <= 1;
          otherwise 8*N.
          For optimal efficiency, LWORK >= (NB+3)*N,
          where NB is the max of the blocksize for SSYTRD and SORMTR
          returned by ILAENV.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (5*N)
[out]IFAIL
          IFAIL is INTEGER array, dimension (N)
          If JOBZ = 'V', then if INFO = 0, the first M elements of
          IFAIL are zero.  If INFO > 0, then IFAIL contains the
          indices of the eigenvectors that failed to converge.
          If JOBZ = 'N', then IFAIL is not referenced.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, then i eigenvectors failed to converge.
                Their indices are stored in array IFAIL.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 250 of file ssyevx.f.

253*
254* -- LAPACK driver routine --
255* -- LAPACK is a software package provided by Univ. of Tennessee, --
256* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
257*
258* .. Scalar Arguments ..
259 CHARACTER JOBZ, RANGE, UPLO
260 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
261 REAL ABSTOL, VL, VU
262* ..
263* .. Array Arguments ..
264 INTEGER IFAIL( * ), IWORK( * )
265 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
266* ..
267*
268* =====================================================================
269*
270* .. Parameters ..
271 REAL ZERO, ONE
272 parameter( zero = 0.0e+0, one = 1.0e+0 )
273* ..
274* .. Local Scalars ..
275 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
276 $ WANTZ
277 CHARACTER ORDER
278 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
279 $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
280 $ ITMP1, J, JJ, LLWORK, LLWRKN, LWKMIN,
281 $ LWKOPT, NB, NSPLIT
282 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
283 $ SIGMA, SMLNUM, TMP1, VLL, VUU
284* ..
285* .. External Functions ..
286 LOGICAL LSAME
287 INTEGER ILAENV
288 REAL SLAMCH, SLANSY, SROUNDUP_LWORK
290* ..
291* .. External Subroutines ..
292 EXTERNAL scopy, slacpy, sorgtr, sormtr, sscal, sstebz,
294* ..
295* .. Intrinsic Functions ..
296 INTRINSIC max, min, sqrt
297* ..
298* .. Executable Statements ..
299*
300* Test the input parameters.
301*
302 lower = lsame( uplo, 'L' )
303 wantz = lsame( jobz, 'V' )
304 alleig = lsame( range, 'A' )
305 valeig = lsame( range, 'V' )
306 indeig = lsame( range, 'I' )
307 lquery = ( lwork.EQ.-1 )
308*
309 info = 0
310 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
311 info = -1
312 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
313 info = -2
314 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
315 info = -3
316 ELSE IF( n.LT.0 ) THEN
317 info = -4
318 ELSE IF( lda.LT.max( 1, n ) ) THEN
319 info = -6
320 ELSE
321 IF( valeig ) THEN
322 IF( n.GT.0 .AND. vu.LE.vl )
323 $ info = -8
324 ELSE IF( indeig ) THEN
325 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
326 info = -9
327 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
328 info = -10
329 END IF
330 END IF
331 END IF
332 IF( info.EQ.0 ) THEN
333 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
334 info = -15
335 END IF
336 END IF
337*
338 IF( info.EQ.0 ) THEN
339 IF( n.LE.1 ) THEN
340 lwkmin = 1
341 work( 1 ) = sroundup_lwork(lwkmin)
342 ELSE
343 lwkmin = 8*n
344 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
345 nb = max( nb, ilaenv( 1, 'SORMTR', uplo, n, -1, -1, -1 ) )
346 lwkopt = max( lwkmin, ( nb + 3 )*n )
347 work( 1 ) = sroundup_lwork(lwkopt)
348 END IF
349*
350 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
351 $ info = -17
352 END IF
353*
354 IF( info.NE.0 ) THEN
355 CALL xerbla( 'SSYEVX', -info )
356 RETURN
357 ELSE IF( lquery ) THEN
358 RETURN
359 END IF
360*
361* Quick return if possible
362*
363 m = 0
364 IF( n.EQ.0 ) THEN
365 RETURN
366 END IF
367*
368 IF( n.EQ.1 ) THEN
369 IF( alleig .OR. indeig ) THEN
370 m = 1
371 w( 1 ) = a( 1, 1 )
372 ELSE
373 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
374 m = 1
375 w( 1 ) = a( 1, 1 )
376 END IF
377 END IF
378 IF( wantz )
379 $ z( 1, 1 ) = one
380 RETURN
381 END IF
382*
383* Get machine constants.
384*
385 safmin = slamch( 'Safe minimum' )
386 eps = slamch( 'Precision' )
387 smlnum = safmin / eps
388 bignum = one / smlnum
389 rmin = sqrt( smlnum )
390 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
391*
392* Scale matrix to allowable range, if necessary.
393*
394 iscale = 0
395 abstll = abstol
396 IF( valeig ) THEN
397 vll = vl
398 vuu = vu
399 END IF
400 anrm = slansy( 'M', uplo, n, a, lda, work )
401 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
402 iscale = 1
403 sigma = rmin / anrm
404 ELSE IF( anrm.GT.rmax ) THEN
405 iscale = 1
406 sigma = rmax / anrm
407 END IF
408 IF( iscale.EQ.1 ) THEN
409 IF( lower ) THEN
410 DO 10 j = 1, n
411 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
412 10 CONTINUE
413 ELSE
414 DO 20 j = 1, n
415 CALL sscal( j, sigma, a( 1, j ), 1 )
416 20 CONTINUE
417 END IF
418 IF( abstol.GT.0 )
419 $ abstll = abstol*sigma
420 IF( valeig ) THEN
421 vll = vl*sigma
422 vuu = vu*sigma
423 END IF
424 END IF
425*
426* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
427*
428 indtau = 1
429 inde = indtau + n
430 indd = inde + n
431 indwrk = indd + n
432 llwork = lwork - indwrk + 1
433 CALL ssytrd( uplo, n, a, lda, work( indd ), work( inde ),
434 $ work( indtau ), work( indwrk ), llwork, iinfo )
435*
436* If all eigenvalues are desired and ABSTOL is less than or equal to
437* zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
438* some eigenvalue, then try SSTEBZ.
439*
440 test = .false.
441 IF( indeig ) THEN
442 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
443 test = .true.
444 END IF
445 END IF
446 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
447 CALL scopy( n, work( indd ), 1, w, 1 )
448 indee = indwrk + 2*n
449 IF( .NOT.wantz ) THEN
450 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
451 CALL ssterf( n, w, work( indee ), info )
452 ELSE
453 CALL slacpy( 'A', n, n, a, lda, z, ldz )
454 CALL sorgtr( uplo, n, z, ldz, work( indtau ),
455 $ work( indwrk ), llwork, iinfo )
456 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
457 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
458 $ work( indwrk ), info )
459 IF( info.EQ.0 ) THEN
460 DO 30 i = 1, n
461 ifail( i ) = 0
462 30 CONTINUE
463 END IF
464 END IF
465 IF( info.EQ.0 ) THEN
466 m = n
467 GO TO 40
468 END IF
469 info = 0
470 END IF
471*
472* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
473*
474 IF( wantz ) THEN
475 order = 'B'
476 ELSE
477 order = 'E'
478 END IF
479 indibl = 1
480 indisp = indibl + n
481 indiwo = indisp + n
482 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
483 $ work( indd ), work( inde ), m, nsplit, w,
484 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
485 $ iwork( indiwo ), info )
486*
487 IF( wantz ) THEN
488 CALL sstein( n, work( indd ), work( inde ), m, w,
489 $ iwork( indibl ), iwork( indisp ), z, ldz,
490 $ work( indwrk ), iwork( indiwo ), ifail, info )
491*
492* Apply orthogonal matrix used in reduction to tridiagonal
493* form to eigenvectors returned by SSTEIN.
494*
495 indwkn = inde
496 llwrkn = lwork - indwkn + 1
497 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
498 $ ldz, work( indwkn ), llwrkn, iinfo )
499 END IF
500*
501* If matrix was scaled, then rescale eigenvalues appropriately.
502*
503 40 CONTINUE
504 IF( iscale.EQ.1 ) THEN
505 IF( info.EQ.0 ) THEN
506 imax = m
507 ELSE
508 imax = info - 1
509 END IF
510 CALL sscal( imax, one / sigma, w, 1 )
511 END IF
512*
513* If eigenvalues are not in order, then sort them, along with
514* eigenvectors.
515*
516 IF( wantz ) THEN
517 DO 60 j = 1, m - 1
518 i = 0
519 tmp1 = w( j )
520 DO 50 jj = j + 1, m
521 IF( w( jj ).LT.tmp1 ) THEN
522 i = jj
523 tmp1 = w( jj )
524 END IF
525 50 CONTINUE
526*
527 IF( i.NE.0 ) THEN
528 itmp1 = iwork( indibl+i-1 )
529 w( i ) = w( j )
530 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
531 w( j ) = tmp1
532 iwork( indibl+j-1 ) = itmp1
533 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
534 IF( info.NE.0 ) THEN
535 itmp1 = ifail( i )
536 ifail( i ) = ifail( j )
537 ifail( j ) = itmp1
538 END IF
539 END IF
540 60 CONTINUE
541 END IF
542*
543* Set WORK(1) to optimal workspace size.
544*
545 work( 1 ) = sroundup_lwork(lwkopt)
546*
547 RETURN
548*
549* End of SSYEVX
550*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ssytrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
SSYTRD
Definition ssytrd.f:192
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
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:273
subroutine sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:174
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
subroutine sorgtr(uplo, n, a, lda, tau, work, lwork, info)
SORGTR
Definition sorgtr.f:123
subroutine sormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
SORMTR
Definition sormtr.f:172
Here is the call graph for this function:
Here is the caller graph for this function: