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

◆ ssyevx_2stage()

subroutine ssyevx_2stage ( 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_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
 SSYEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric matrix A using the 2stage technique for
 the reduction to tridiagonal.  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.
                  Not available in this release.
[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
          If JOBZ = 'N' and N > 1, LWORK must be queried.
                                   LWORK = MAX(1, 8*N, dimension) where
                                   dimension = max(stage1,stage2) + (KD+1)*N + 3*N
                                             = N*KD + N*max(KD+1,FACTOPTNB)
                                               + max(2*KD*KD, KD*NTHREADS)
                                               + (KD+1)*N + 3*N
                                   where KD is the blocking size of the reduction,
                                   FACTOPTNB is the blocking used by the QR or LQ
                                   algorithm, usually FACTOPTNB=128 is a good choice
                                   NTHREADS is the number of threads used when
                                   openMP compilation is enabled, otherwise =1.
          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

          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.
Further Details:
  All details about the 2stage techniques are available in:

  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  Parallel reduction to condensed forms for symmetric eigenvalue problems
  using aggregated fine-grained and memory-aware kernels. In Proceedings
  of 2011 International Conference for High Performance Computing,
  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  Article 8 , 11 pages.
  http://doi.acm.org/10.1145/2063384.2063394

  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  An improved parallel singular value algorithm and its implementation
  for multicore hardware, In Proceedings of 2013 International Conference
  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  Denver, Colorado, USA, 2013.
  Article 90, 12 pages.
  http://doi.acm.org/10.1145/2503210.2503292

  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  A novel hybrid CPU-GPU generalized eigensolver for electronic structure
  calculations based on fine-grained memory aware tasks.
  International Journal of High Performance Computing Applications.
  Volume 28 Issue 2, Pages 196-209, May 2014.
  http://hpc.sagepub.com/content/28/2/196

Definition at line 297 of file ssyevx_2stage.f.

300*
301 IMPLICIT NONE
302*
303* -- LAPACK driver routine --
304* -- LAPACK is a software package provided by Univ. of Tennessee, --
305* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306*
307* .. Scalar Arguments ..
308 CHARACTER JOBZ, RANGE, UPLO
309 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
310 REAL ABSTOL, VL, VU
311* ..
312* .. Array Arguments ..
313 INTEGER IFAIL( * ), IWORK( * )
314 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
315* ..
316*
317* =====================================================================
318*
319* .. Parameters ..
320 REAL ZERO, ONE
321 parameter( zero = 0.0e+0, one = 1.0e+0 )
322* ..
323* .. Local Scalars ..
324 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
325 $ WANTZ
326 CHARACTER ORDER
327 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
328 $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
329 $ ITMP1, J, JJ, LLWORK, LLWRKN,
330 $ NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
331 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
332 $ SIGMA, SMLNUM, TMP1, VLL, VUU
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 INTEGER ILAENV2STAGE
337 REAL SLAMCH, SLANSY, SROUNDUP_LWORK
338 EXTERNAL lsame, slamch, slansy, ilaenv2stage,
340* ..
341* .. External Subroutines ..
342 EXTERNAL scopy, slacpy, sorgtr, sormtr, sscal, sstebz,
345* ..
346* .. Intrinsic Functions ..
347 INTRINSIC max, min, sqrt
348* ..
349* .. Executable Statements ..
350*
351* Test the input parameters.
352*
353 lower = lsame( uplo, 'L' )
354 wantz = lsame( jobz, 'V' )
355 alleig = lsame( range, 'A' )
356 valeig = lsame( range, 'V' )
357 indeig = lsame( range, 'I' )
358 lquery = ( lwork.EQ.-1 )
359*
360 info = 0
361 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
362 info = -1
363 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
364 info = -2
365 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
366 info = -3
367 ELSE IF( n.LT.0 ) THEN
368 info = -4
369 ELSE IF( lda.LT.max( 1, n ) ) THEN
370 info = -6
371 ELSE
372 IF( valeig ) THEN
373 IF( n.GT.0 .AND. vu.LE.vl )
374 $ info = -8
375 ELSE IF( indeig ) THEN
376 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
377 info = -9
378 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
379 info = -10
380 END IF
381 END IF
382 END IF
383 IF( info.EQ.0 ) THEN
384 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
385 info = -15
386 END IF
387 END IF
388*
389 IF( info.EQ.0 ) THEN
390 IF( n.LE.1 ) THEN
391 lwmin = 1
392 work( 1 ) = sroundup_lwork(lwmin)
393 ELSE
394 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz,
395 $ n, -1, -1, -1 )
396 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz,
397 $ n, kd, -1, -1 )
398 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz,
399 $ n, kd, ib, -1 )
400 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz,
401 $ n, kd, ib, -1 )
402 lwmin = max( 8*n, 3*n + lhtrd + lwtrd )
403 work( 1 ) = lwmin
404 END IF
405*
406 IF( lwork.LT.lwmin .AND. .NOT.lquery )
407 $ info = -17
408 END IF
409*
410 IF( info.NE.0 ) THEN
411 CALL xerbla( 'SSYEVX_2STAGE', -info )
412 RETURN
413 ELSE IF( lquery ) THEN
414 RETURN
415 END IF
416*
417* Quick return if possible
418*
419 m = 0
420 IF( n.EQ.0 ) THEN
421 RETURN
422 END IF
423*
424 IF( n.EQ.1 ) THEN
425 IF( alleig .OR. indeig ) THEN
426 m = 1
427 w( 1 ) = a( 1, 1 )
428 ELSE
429 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
430 m = 1
431 w( 1 ) = a( 1, 1 )
432 END IF
433 END IF
434 IF( wantz )
435 $ z( 1, 1 ) = one
436 RETURN
437 END IF
438*
439* Get machine constants.
440*
441 safmin = slamch( 'Safe minimum' )
442 eps = slamch( 'Precision' )
443 smlnum = safmin / eps
444 bignum = one / smlnum
445 rmin = sqrt( smlnum )
446 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
447*
448* Scale matrix to allowable range, if necessary.
449*
450 iscale = 0
451 abstll = abstol
452 IF( valeig ) THEN
453 vll = vl
454 vuu = vu
455 END IF
456 anrm = slansy( 'M', uplo, n, a, lda, work )
457 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
458 iscale = 1
459 sigma = rmin / anrm
460 ELSE IF( anrm.GT.rmax ) THEN
461 iscale = 1
462 sigma = rmax / anrm
463 END IF
464 IF( iscale.EQ.1 ) THEN
465 IF( lower ) THEN
466 DO 10 j = 1, n
467 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
468 10 CONTINUE
469 ELSE
470 DO 20 j = 1, n
471 CALL sscal( j, sigma, a( 1, j ), 1 )
472 20 CONTINUE
473 END IF
474 IF( abstol.GT.0 )
475 $ abstll = abstol*sigma
476 IF( valeig ) THEN
477 vll = vl*sigma
478 vuu = vu*sigma
479 END IF
480 END IF
481*
482* Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
483*
484 indtau = 1
485 inde = indtau + n
486 indd = inde + n
487 indhous = indd + n
488 indwrk = indhous + lhtrd
489 llwork = lwork - indwrk + 1
490*
491 CALL ssytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
492 $ work( inde ), work( indtau ), work( indhous ),
493 $ lhtrd, work( indwrk ), llwork, iinfo )
494*
495* If all eigenvalues are desired and ABSTOL is less than or equal to
496* zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
497* some eigenvalue, then try SSTEBZ.
498*
499 test = .false.
500 IF( indeig ) THEN
501 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
502 test = .true.
503 END IF
504 END IF
505 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
506 CALL scopy( n, work( indd ), 1, w, 1 )
507 indee = indwrk + 2*n
508 IF( .NOT.wantz ) THEN
509 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
510 CALL ssterf( n, w, work( indee ), info )
511 ELSE
512 CALL slacpy( 'A', n, n, a, lda, z, ldz )
513 CALL sorgtr( uplo, n, z, ldz, work( indtau ),
514 $ work( indwrk ), llwork, iinfo )
515 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
516 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
517 $ work( indwrk ), info )
518 IF( info.EQ.0 ) THEN
519 DO 30 i = 1, n
520 ifail( i ) = 0
521 30 CONTINUE
522 END IF
523 END IF
524 IF( info.EQ.0 ) THEN
525 m = n
526 GO TO 40
527 END IF
528 info = 0
529 END IF
530*
531* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
532*
533 IF( wantz ) THEN
534 order = 'B'
535 ELSE
536 order = 'E'
537 END IF
538 indibl = 1
539 indisp = indibl + n
540 indiwo = indisp + n
541 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
542 $ work( indd ), work( inde ), m, nsplit, w,
543 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
544 $ iwork( indiwo ), info )
545*
546 IF( wantz ) THEN
547 CALL sstein( n, work( indd ), work( inde ), m, w,
548 $ iwork( indibl ), iwork( indisp ), z, ldz,
549 $ work( indwrk ), iwork( indiwo ), ifail, info )
550*
551* Apply orthogonal matrix used in reduction to tridiagonal
552* form to eigenvectors returned by SSTEIN.
553*
554 indwkn = inde
555 llwrkn = lwork - indwkn + 1
556 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
557 $ ldz, work( indwkn ), llwrkn, iinfo )
558 END IF
559*
560* If matrix was scaled, then rescale eigenvalues appropriately.
561*
562 40 CONTINUE
563 IF( iscale.EQ.1 ) THEN
564 IF( info.EQ.0 ) THEN
565 imax = m
566 ELSE
567 imax = info - 1
568 END IF
569 CALL sscal( imax, one / sigma, w, 1 )
570 END IF
571*
572* If eigenvalues are not in order, then sort them, along with
573* eigenvectors.
574*
575 IF( wantz ) THEN
576 DO 60 j = 1, m - 1
577 i = 0
578 tmp1 = w( j )
579 DO 50 jj = j + 1, m
580 IF( w( jj ).LT.tmp1 ) THEN
581 i = jj
582 tmp1 = w( jj )
583 END IF
584 50 CONTINUE
585*
586 IF( i.NE.0 ) THEN
587 itmp1 = iwork( indibl+i-1 )
588 w( i ) = w( j )
589 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
590 w( j ) = tmp1
591 iwork( indibl+j-1 ) = itmp1
592 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
593 IF( info.NE.0 ) THEN
594 itmp1 = ifail( i )
595 ifail( i ) = ifail( j )
596 ifail( j ) = itmp1
597 END IF
598 END IF
599 60 CONTINUE
600 END IF
601*
602* Set WORK(1) to optimal workspace size.
603*
604 work( 1 ) = sroundup_lwork(lwmin)
605*
606 RETURN
607*
608* End of SSYEVX_2STAGE
609*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ssytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
SSYTRD_2STAGE
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
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: